{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "56eede84",
   "metadata": {},
   "source": [
    "# Solving nonlinear diffusion-reaction system with different priors\n",
    "\n",
    "We consider a one-dimensional nonlinear diffusion-reaction equation, with unknown constant $k_r$:\n",
    "\n",
    "$$D\\partial_x^2 u - k_r u^3 = f, x\\in[-1, 1], D=0.01.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "389d971b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import neuraluq as neuq\n",
    "import neuraluq.variables as neuq_vars\n",
    "from neuraluq.config import tf\n",
    "\n",
    "import numpy as np\n",
    "import scipy.io as sio\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50d38f55",
   "metadata": {},
   "source": [
    "#### Common functions for loading data, visualization, building Bayesian models and PDE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ea30ef6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot1d(x, y, x_test, y_test, y_samples, xlim=None, ylim=None, title=\"\", ylabel=\"\"):\n",
    "    y_mean = np.mean(y_samples, axis=0)\n",
    "    y_std = np.std(y_samples, axis=0)\n",
    "\n",
    "    # plt.figure(figsize=(8, 6), dpi=80)\n",
    "    plt.plot(x, y, \"k.\", markersize=10, label=\"Training data\")\n",
    "    plt.plot(x_test, y_test, \"k-\", label=\"Exact\")\n",
    "    plt.plot(x_test, y_mean, \"r--\", label=\"Mean\")\n",
    "    plt.fill_between(\n",
    "        x_test.ravel(),\n",
    "        y_mean + 2 * y_std,\n",
    "        y_mean - 2 * y_std,\n",
    "        alpha=0.3,\n",
    "        facecolor=\"c\",\n",
    "        label=\"2 stds\",\n",
    "    )\n",
    "    plt.legend()\n",
    "    plt.xlabel(\"$x$\")\n",
    "    plt.ylabel(ylabel)\n",
    "    plt.xlim(xlim)\n",
    "    plt.ylim(ylim)\n",
    "    plt.title(title)\n",
    "    plt.show()\n",
    "\n",
    "def load_data(noise):\n",
    "    data = sio.loadmat(\"../../dataset/1D_different_priors.mat\")\n",
    "    x_test, u_test, f_test = data[\"x_test\"], data[\"u_test\"], data[\"f_test\"]\n",
    "    idx = np.random.choice(100, 5, replace=False)\n",
    "    x_u_train, u_train = x_test[idx, ...], u_test[idx, ...]\n",
    "    idx = np.arange(100)[::6]\n",
    "    x_f_train, f_train = x_test[idx, ...], f_test[idx, ...]\n",
    "    f_train = f_train + noise*np.random.normal(size=f_train.shape)\n",
    "    u_train = u_train + noise*np.random.normal(size=u_train.shape)\n",
    "    return x_f_train, f_train, x_u_train, u_train, x_test, u_test, f_test\n",
    "\n",
    "\n",
    "def pde_fn(x, u, k):\n",
    "    D = 0.01\n",
    "    u_x = tf.gradients(u, x)[0]\n",
    "    u_xx = tf.gradients(u_x, x)[0]\n",
    "    return D * u_xx - k * u**3\n",
    "\n",
    "\n",
    "def build_model(process_u, process_k, burnin_steps=1000):\n",
    "    likelihood_f = neuq.likelihoods.Normal(\n",
    "        inputs=x_f_train,\n",
    "        targets=f_train,\n",
    "        processes=[process_u, process_k],\n",
    "        pde=pde_fn,\n",
    "        sigma=noise,\n",
    "    )\n",
    "    likelihood_u = neuq.likelihoods.Normal(\n",
    "        inputs=x_u_train,\n",
    "        targets=u_train,\n",
    "        processes=[process_u],\n",
    "        sigma=noise,\n",
    "    )\n",
    "\n",
    "    model = neuq.models.Model(\n",
    "        processes=[process_u, process_k],\n",
    "        likelihoods=[likelihood_u, likelihood_f],\n",
    "    )\n",
    "\n",
    "    method = neuq.inferences.HMC(\n",
    "        num_samples=1000,\n",
    "        num_burnin=burnin_steps,\n",
    "        init_time_step=0.01,\n",
    "        leapfrog_step=50,\n",
    "        seed=6666,\n",
    "    )\n",
    "    model.compile(method=method)\n",
    "    samples, results = model.run()\n",
    "    print(\"Acceptance rate: \", np.mean(results))\n",
    "    return samples, model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "358a5ac2",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(666)\n",
    "noise = 0.1\n",
    "x_f_train, f_train, x_u_train, u_train, x_test, u_test, f_test = load_data(noise)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e563cf08",
   "metadata": {},
   "source": [
    "#### Different priors for the same surrogate\n",
    "\n",
    "In this problem, we assume $k_r$ is unknown and constant, and we give it different priors, Normal, Half Normal and Log Normal:\n",
    "\n",
    "1. $k_r\\sim N(0, 1^2)$\n",
    "2. $k_r \\sim HalfN(0, 1^2)$\n",
    "3. $k_r \\sim LogN(0, 1^2)$\n",
    "\n",
    "Normal for BNN is a built-in distribution, and hence we need to customize the other two.\n",
    "\n",
    "The exact value of $k_r$ used to generator data is $0.2$.\n",
    "\n",
    "Choosing an appropriate prior is crucial to the estimation, for some times. For example, if the prior knowledge tells us that a quantity is non-negative, then Normal as the prior distribution is, most of times, inferior to some other distributions, like Half Normal. For strictly positive quantities, like standard deviation, Gamma and Log Normal distributions may be better choices.\n",
    "\n",
    "For this example, if we know that $k_r$ is non-negative or strictly positive, choosing Half Normal or Log Normal, together with feasible initial states for Markov chain, will guarantee that the posterior samples are physically reasonable."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "994608f1",
   "metadata": {},
   "source": [
    "#### Normal prior for $k_r$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a8ebcd2f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Supporting backend tensorflow.compat.v1\n",
      "\n",
      "Compiling a MCMC method\n",
      "\n",
      "sampling from posterior distribution ...\n",
      "\n",
      "Finished sampling from posterior distribution ...\n",
      "\n",
      "Acceptance rate:  0.614\n"
     ]
    }
   ],
   "source": [
    "layers = [1, 50, 50, 50, 1]\n",
    "prior_u = neuq.variables.fnn.Samplable(layers=layers, mean=0, sigma=1) # Normal\n",
    "process_u = neuq.process.Process(\n",
    "    surrogate=neuq.surrogates.FNN(layers=layers),\n",
    "    prior=prior_u,\n",
    ")\n",
    "prior_k = neuq.variables.const.Samplable(mean=0, sigma=1)\n",
    "\n",
    "process_k = neuq.process.Process(\n",
    "    surrogate=neuq.surrogates.Identity(),\n",
    "    prior=prior_k,\n",
    ")\n",
    "# use burnin steps to tune to good acceptance rate, with random seed fixed\n",
    "samples, model = build_model(process_u, process_k, burnin_steps=800)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3fd3ec48",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEaCAYAAAAPGBBTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABTVUlEQVR4nO3dd3hUVfrA8e+ZmUySmVR6DV1AWhIjSJEiFppipSgriC6KsJa1rB3BvohiV+xlRUBF8SeIoLAoRQi9Q+ghdNLbTGbO748M2QHSMy3h/TxPnmTu3Ln35GYy7z3tPUprjRBCCFESg78LIIQQIrBJoBBCCFEqCRRCCCFKJYFCCCFEqSRQCCGEKJUECiGEEKWSQCGEEKJUEiiECHBKqa5KqZVKqWVKqZlKqSB/l0lcWCRQCBH4DgFXaK17A/uBof4tjrjQSKAQogqUUi8ppR7w5jm01ke01rmuhzbA6Tr3aqVUB2+eWwiQQCFEiZRS0UoprZRaec7295VSryul6gK3Ax+4PVdLKTVXKZWtlDqglLq1nOeaqJRKVErlK6U+K2GfZsDVwE+uTa8CUyrxqwlRIRIohChZLHAUuFgp1cBtexywARgDzHe72wd4h8K7/vrAbcB75bzrTwGeBz4p7kmlVATwJTBGa213bZ4H9DunbEJ4nAQKIUoWCyQCi3D1CyiljEAnYD0wEPjvmZ2VUlbgJuBprXWW1vpPCj/M/1bWibTW32utfwBOnfucUsoEfANM1lrvdHtNHrAWuKZyv54Q5SOBQoiSnak5/ABc79rWjsL/m+0UBoydbvtfBBRorXe5bdsIVLUfYSTQDXhaKbVUKTXc7bntQJcqHl+IUpn8XQAhAlgs8CPwO/C+UirctW2r1tqulIoCMt32DwMyzjlGOhBelUJorb+ksNmpOJlAw6ocX4iySI1CiGIopYKB9sAGrXUqsJrCpqYztQyAVM4OAllAxDmHiuDsYOJp4UCaF48vhAQKIUrQEcgB9roe/0Bh81Mchf0TAJsobG46YxdgUkq1cdvWBdjqxXK2p7B5SwivkUAhRPHigE36f0tAzgMGcXaNYj7Q58wLtNbZwPfAFKWUVSnVk8JO8KJmI6XUZ8UNf1VKmZRSIYARMCqlQlyd2CVy7X8JhZ3tQniNBAohihfL/wICWuv9FM6KjuJ/d/BfAIOUUqFur7sXCAWOAzOB8Vpr9xpFU2B5Med7CsgFHgNGuX5+qowyXgss1VqnlP3rCFF5StbMFqLylFIvAse11tPLsa+ZwiDT2W0uRFXO/Rdwp9Z6S1WPJURpJFAIIYQolTQ9CSGEKJXfAoVSqqlSaolSaptSaqtS6v5i9lFKqTeVUklKqU1KqXh/lFUIIS5k/pxwVwA8pLVe55rItFYptUhrvc1tn4FAG9dXN+A913chhBA+4rcahSt18jrXz5kUpiJofM5uQ4EvdKFVQJRSSmahCiGEDwVECg+lVHMKx6f/dc5TjSlctOWMZNe2I8UcYxwwDsBqtV7Srl07r5RVCCFqorVr157UWtct7jm/BwqlVBjwHfCA1vrcPDnlprWeAcwASEhI0ImJiR4qoRBC1HxKqQMlPefXUU+utX+/A/6jtf6+mF0OUzhB6Ywmrm1CCCF8xJ+jnhTwMbBda/1aCbvNA253jX66DEjXWp/X7CSEEMJ7/Nn01JPCBV02K6U2uLY9AcQAaK3fpzCXziAgicIEbXf4vphCCHFh81ugcK3+pcrYRwMTfFMiIYQ32O12kpOTycvL83dRBBASEkKTJk0ICgoq92v83pkthKjZkpOTCQ8Pp3nz5hS2OAt/0Vpz6tQpkpOTadGiRblfJyk8hBBelZeXR+3atSVIBAClFLVr165w7U4ChRDC6yRIBI7K/C0kUAghhCiVBAohRI126tQpYmNjiY2NpUGDBjRu3Ljosc1mK/W1iYmJ3HfffWWeo0ePHp4q7ln69u1LWZOHp0+fTk5OjlfOf4Z0ZgshAorD4WDBggWsX7+euLg4Bg4ciNForPTxateuzYYNGwB49tlnCQsL4+GHHy56vqCgAJOp+I/ChIQEEhISyjzHihUrKl2+qpo+fTqjRo3CYrF47RxSoxBCBAyHw8E111zDyJEjmTRpEiNHjuSaa67B4XB49DxjxozhnnvuoVu3bjz66KOsXr2a7t27ExcXR48ePdi5cycAS5cuZciQIUBhkBk7dix9+/alZcuWvPnmm0XHCwsLK9q/b9++3HzzzbRr147bbruNM4vDzZ8/n3bt2nHJJZdw3333FR3XXW5uLiNGjKB9+/bccMMN5ObmFj03fvx4EhIS6NChA5MmTQLgzTffJCUlhX79+tGvX78S96sqqVEIIQLGggUL+Ouvv8jKygIgKyuLv/76iwULFhT7wVoVycnJrFixAqPRSEZGBn/88Qcmk4nFixfzxBNP8N133533mh07drBkyRIyMzNp27Yt48ePP28+wvr169m6dSuNGjWiZ8+eLF++nISEBO6++26WLVtGixYtGDlyZLFleu+997BYLGzfvp1NmzYRH/+/JXheeOEFatWqhcPhoH///mzatIn77ruP1157jSVLllCnTp0S9+vcuXOVrpXUKIQQAWP9+vVkZ2eftS07O7uo6ciTbrnllqImrfT0dG655RY6duzIgw8+yNatW4t9zeDBgwkODqZOnTrUq1ePY8eOnbdP165dadKkCQaDgdjYWPbv38+OHTto2bJl0dyFkgLFsmXLGDVqFACdO3c+6wN+9uzZxMfHExcXx9atW9m2bVuxxyjvfhUhgUIIETDi4uKwWq1nbbNarcTGxnr8XO7nefrpp+nXrx9btmzhp59+KnGeQXBwcNHPRqORgoKCSu1TUfv27ePVV1/lt99+Y9OmTQwePLjYMpZ3v4qSQCGECBgDBw6kW7duhIWFoZQiLCyMbt26MXDgQK+eNz09ncaNC9dN++yzzzx+/LZt27J37172798PwKxZs4rdr3fv3nz99dcAbNmyhU2bNgGQkZGB1WolMjKSY8eOsWDBgqLXhIeHk5mZWeZ+VSF9FEKIgGE0Glm4cCELFixgw4YNxMbGVnnUU3k8+uijjB49mueff57Bgwd7/PihoaG8++67DBgwAKvVyqWXXlrsfuPHj+eOO+6gffv2tG/fnksuuQSALl26EBcXR7t27WjatCk9e/Yses24ceMYMGAAjRo1YsmSJSXuVxXqTI98TSILFwkROLZv30779u39XQy/y8rKIiwsDK01EyZMoE2bNjz44IN+KUtxfxOl1FqtdbFjgaXpSQghfODDDz8kNjaWDh06kJ6ezt133+3vIpWbND0JIYQPPPjgg36rQVSV1CiEEEKUSgKFEEKIUkmgEEIIUSq/Bgql1CdKqeNKqS0lPN9XKZWulNrg+nrGm+WxOZ3ePLwQQlRL/q5RfAYMKGOfP7TWsa6vKd4sTI7DQZrd7s1TCCH8wGg0FqUWj42N5eWXX/bYsTds2MD8+fM9drxA5NdRT1rrZUqp5v4sw7kO22xEVWDRcSFE4AsNDfVKvigoDBSJiYkMGjTIK8cPBP6uUZRHd6XURqXUAqVUB2+fLCU/n5o4CVEIcbb09HTatm1blFJ85MiRfPjhh0DJqbrXrFlDjx496NKlC127diU9PZ1nnnmGWbNmERsbW2Jqjuou0OdRrAOaaa2zlFKDgB+ANsXtqJQaB4wDiImJqfQJ85xOTtnt1DGbK30MIUTxHnjgAY/f2cfGxjJ9+vRS98nNzT0rseDjjz/O8OHDefvttxkzZgz3338/qamp/P3vfweKT9Xdrl07hg8fzqxZs7j00kvJyMjAYrEwZcoUEhMTefvttz36ewWSgA4UWusMt5/nK6XeVUrV0VqfLGbfGcAMKEzhUZXzHrbZJFAIUYOU1PR01VVXMWfOHCZMmMDGjRuLts+ePZsZM2ZQUFDAkSNH2LZtG0opGjZsWJSnKSIiwlfF97uADhRKqQbAMa21Vkp1pbCp7JS3z3skP59OVisGpbx9KiEuKGXd+fua0+lk+/btWCwWUlNTadKkSVGq7jVr1hAdHc2YMWM8kqq7OvP38NiZwEqgrVIqWSl1p1LqHqXUPa5dbga2KKU2Am8CI7QPOhDsWnO8jEXXhRDV3+uvv0779u35+uuvueOOO7Db7SWm6m7bti1HjhxhzZo1AGRmZlJQUHBWmu+ayt+jnopf5ul/z78N+KXhL8Vmo4HbAiRCiOrr3D6KAQMGcMcdd/DRRx+xevVqwsPD6d27N88//zyTJ08uNlW32Wxm1qxZ/OMf/yA3N5fQ0FAWL15Mv379ePnll4mNjS3q+6hpJM24mzS7nT/S0wEwKsVV0dEEGarDwDAhApekGQ88kmbcQxxac0San4QQQgJFaQ5e4B1YQggBEihKlVpQQJYHFkYXQojqTAJFGQ7l5/u7CEII4VcSKMpwSFJ6CCEucBIoypDvdHJcMsoKIS5gEijK4ZB0agtRrSmlGDVqVNHjgoIC6taty5AhQ/xYqupDAkU5HLPbyZdFjYSotqxWK1u2bCE3NxeARYsW0bhxYz+XqvqQQFEOTq1lqKwQ1dygQYP4+eefAZg5cyYjR/4vMUR2djZjx46la9euxMXF8eOPPwKwf/9+Lr/8cuLj44mPj2fFihUALF26lL59+3LzzTfTrl07brvtthrdlxnQSQEDyYG8PFqHhqIkUaAQVdO37/nbhg2De++FnBwobgGgMWMKv06ehJtvPvu5pUvLddoRI0YwZcoUhgwZwqZNmxg7dix//PEHUJhW/IorruCTTz4hLS2Nrl27cuWVV1KvXj0WLVpESEgIu3fvZuTIkZzJ+rB+/Xq2bt1Ko0aN6NmzJ8uXL6dXr17lvgzViQSKcsp1Ojlht1NP0o8LUS117tyZ/fv3M3PmzPNWo/v111+ZN28er776KgB5eXkcPHiQRo0aMXHiRDZs2IDRaGTXrl1Fr+natStNmjQBCtfE2L9/vwQKAfvz8iRQCFFVpdUALJbSn69Tp9w1iOJcd911PPzwwyxdupRTp/63YoHWmu+++462bduetf+zzz5L/fr12bhxI06nk5CQkKLngt2ShhqNRgpq8ORc6aOogGM2G7kOh7+LIYSopLFjxzJp0iQ6dep01vZrrrmGt956q6ifYf369UDhcqkNGzbEYDDw5Zdf4rhA//8lUFTQAenUFqLaatKkCffdd995259++mnsdjudO3emQ4cOPP300wDce++9fP7553Tp0oUdO3ZgtVp9XeSAIGnG3binGS9JsMHAldHRsvqdEOUkacYDj6QZ97J8p5MUyf8khLiASKCohD3S/CSEuID4e83sT5RSx5VSW0p4Ximl3lRKJSmlNiml4n1dxuJkFBRwQhY1EkJcIPxdo/gMGFDK8wOBNq6vccB7PihTuexxpQIQQoiazq/zKLTWy5RSzUvZZSjwhS7scV+llIpSSjXUWh/xTQlLdsJuJ6OggAiTfy5hnsPBSbuddIeD9IICFBBhMhFhNNLAbJa1voUQHhPoE+4aA4fcHie7tvk9UEBhrSIuPNyn5zxtt7MvL48j+fmcO17tpCsdutlgoJ3FQkxwsKQcEUJUWaAHinJTSo2jsHmKmJgYn5zzcH4+bS0WLEajV8+jteaYzcbu3FzSyjH70+Z0sikri/15ecSHhRHup1qPEMX56eRJjx7v2jp1Sn3+0KFD3H777Rw7dgylFOPGjeP+++8v17E3bNhASkrKeSk/zmjevDmJiYnUKaMM1V2gt08cBpq6PW7i2nYerfUMrXWC1jqhbt26PimcBrbn5Hjt+E6tSc7L479paazJzCxXkHCXUVDA8vR0TkrHu7iAmUwmpk2bxrZt21i1ahXvvPMO27ZtK9drN2zYwPz5871cwsAX6IFiHnC7a/TTZUB6IPRPuEvJz+e0h1fAy3c62ZmTw+LUVNZnZZFZhbQBdq1ZlZEhiy+JC1bDhg2Jjy8cMBkeHk779u05fPj8+805c+bQsWNHunTpQu/evbHZbDzzzDPMmjWL2NhYZs2axalTp7j66qvp0KEDd911V1HKj+zsbAYPHkyXLl3o2LEjs2bN8unv6G1+bZNQSs0E+gJ1lFLJwCQgCEBr/T4wHxgEJAE5wB3+KWnptmZn0ysyssr9ARkFBezNzeWwzYbTgzPmNbAhKwub1rQKDfXYcYWobvbv38/69evp1q3bec9NmTKFhQsX0rhxY9LS0jCbzUyZMoXExETefvttAO677z569erFM888w88//8zHH38MwC+//EKjRo2K1rtILyPDQ3Xj71FPI8t4XgMTfFScSksrKOBwfj5N3DJLVkRWQQFbc3I47uUmom3Z2Ti1po3F4tXzCBGIsrKyuOmmm5g+fToRERHnPd+zZ0/GjBnDsGHDuPHGG4s9xrJly/j+++8BGDx4MNHR0QB06tSJhx56iH/9618MGTKEyy+/3Hu/iB8EetNTtbE9JwdHBWsBBU4n27KzWZqW5vUgccaOnBx2ebFfRYhAZLfbuemmm7jttttKDALvv/8+zz//PIcOHeKSSy45Kw15WS666CLWrVtHp06deOqpp5gyZYqnih4QJFB4SJ7TydrMzHI3GZ2w2ViSlsae3Nzzhrl6286cHHZLsBAXCK01d955J+3bt+ef//xnifvt2bOHbt26MWXKFOrWrcuhQ4cIDw8nMzOzaJ/evXvz9ddfA7BgwQJSU1MBSElJwWKxMGrUKB555BHWrVvn3V+qBN5K8irjJj3omM3GusxMLgkPL7G/wqE1W7Oz/Z6ufEdODkFK0Vz6LISPlTWc1dOWL1/Ol19+SadOnYiNjQXgxRdfPG/I6yOPPMLu3bvRWtO/f3+6dOlCTEwML7/8MrGxsTz++ONMmjSJkSNH0qFDB3r06FE0FH/z5s088sgjGAwGgoKCeO893yeRKHA6UUrhjcH6kmbcTXnSjJdH4+Bg4sLCzgoWTq1Jzs9nd24uOQG0+El8eDiN3VbqEsLTJM24b2Q7HIQYDBjLMaimomnGpUbhBYfz80krKCDCaCTCZEJRuIxqntPp76KdZ31mJmalqCtLvApRbTm1xu50EuKl1D0SKLwk2+Eg2+HgSIBPdtPA2sxMekdFeX2GuRDCO/K9fBMqndkCu9asycys8KgtIcqrJjZxBwqtNbYKXN/K/C0kUAigcLLfxqwsfxdD1EAhISGcOnVKgoWX2LQu97XVWnPq1ClCKjjnS5qeRJHD+fnUDgqiWSUnDgpRnCZNmpCcnMyJEyf8XZQaKc/pLBqWH2IwYCijMzskJIQmTZpU6BwSKMRZtmZnUzcoSPorhMcEBQXRokULfxejRjpus/FXRkbR435RUYR5IVu0ND2Jszi0ZmNWljQTCFENJPlopU2pUYjznLTbOZCXJ5PxRLmcsts5YbMVrbYIEGowEGowUM9spnFwcLnG9ouKOW23c8rDmatLIoFCFGtbTg71zGZpghLF0lqTYrOxJze3KDi4y3c6SQOO2Gxsy84mJiSEVqGhBMsSvR6zx0e1CZBAIUrg0JrN2dl0KybLpriwZRQUsD4ri4xyLqRl15o9ubkczMujncVCs5AQWaK3ijILCjjqwzlaEt5FiY7bbBzNz/d3MUSA0FqTlJPDH+np5Q4S7uyum4/l6elkVuL14n981TdxhgQKUaot2dkyEU/g0JrVmZlsz8mp8qJaqQUFLEtPJyknRwZNVEKOw8FhH9/ASaAQpcp1OmX9igucQ2vWZGR4dM0Up9Zsz8lhRUYG2QGUJLM62JmT4/OlCSRQiDLtzcsjS5oKLkgOrVmdkcEJL42uOW23szQtjd0eqKlcCDILCkj2Q3OwBApRJqfWbMrO9ncxhB+szczkpJeHYDq1Zoer7yPVR8M9q6udfqrd+zVQKKUGKKV2KqWSlFKPFfP8GKXUCaXUBtfXXf4opygcK3/Qz4stCd/al5vLMR+OrMkoKODP9HTWZWaSK81R50kvKPBbNmq/DY9VShmBd4CrgGRgjVJqntZ62zm7ztJaT/R5AcV5tmVnU99slrHwF4CMggK2+enu9XB+PkdsNlq65l6Y5f0GFK5K6S/+/At0BZK01nu11jbgG2CoH8sjymB3LeMqajaH1hVa/90bnFqTlJvLb6mp7MzJwR6Ai3750kmbzaODCSrKn4GiMXDI7XGya9u5blJKbVJKfauUalrSwZRS45RSiUqpRMlS6T2H8/P9+oYV3rc9O5usAGn6KdCaXTk5LE5NZVt2NnkBUi5fcrrmn/hToNfpfgKaa607A4uAz0vaUWs9Q2udoLVOqFu3rs8KeCHamJXl9RW1hH9kOxzsD8C+qALX7O7Fqamsy8zkhM12wczB2Jub6/fA7c8UHocB9xpCE9e2IlrrU24PPwL+7YNyiTLkOZ2sz8ykW0REwKdisDud5Dqd2LWmQGtqmUwESZt3ibZnZ/t8jH5FaAprtYfz8wkxGGhoNhMdFESk0YjVaDzr/ejUmlynkxyHgxyns2h5YpvrvVCgNUbAbDAQbDAQbjQSbTIRbTJhCpD3SK7DwS4fz8Iujj8DxRqgjVKqBYUBYgRwq/sOSqmGWusjrofXAdt9W0RRkhN2O7tzc7nIYvF3UYqVXlDA3txcUmy2s9rajUrRyGymeUgIUUFBfixh4Dlttwf8Gu/u8pxO9uXlsc9VAzIohXsKS3t5axzn3K0roFZQEI2Dg2lkNvv1xmJrgGRG8Fug0FoXKKUmAgsBI/CJ1nqrUmoKkKi1ngfcp5S6DigATgNj/FVecb6dOTnUMpmoYzb7uyhFTtnt7MjJ4XQJ4/EdWnMoP59D+fm0Dg2lncUS8LUiX9lWzQcqOLXGEw2imsL30Sm7nS1K0dBspnVoKBFeWBCoNEddo78CgV+zx2qt5wPzz9n2jNvPjwOP+7pcovwSMzPpERnp83+ic6W5AkRFZhAnuVJkXxIefsE3Rx3JzydVZt+fx6l1UVNXfbOZthYLkT54r+c4HGwIoDXsJc14DZOdkUHy7t2cPnqU1GPHUAYDTdu0oWnbtkTWqePx89m1ZlVGBj0jI7H6Ye2K4zYbSbm5lV7A5YTdzrL0dLpHRFzQa2/4a8ZvdXLMZuOYzUazkBDaWSxem9/hdA1PLnfTmQ9IoKgB8nNyWL1wIf/9/nvW/fYbBSV8aDZo3pxrbr+dq267jYhatTx3fqeTlenp9IyMJNRHH7ZH8/PZmZtbZrrrArud00ePcjIlhbTjx2l28cU0btXqrH1yHA5WpKfT3U/Bzt+O2WxkXoDDTivrQF4eKfn5tLdaaRYS4vHjb8/JIS3AaneqJg4xS0hI0ImJiRV+XZrdzh/p6V4okXc4CgpY+MUXzJw6lfSTJ6nVoAEDBwygV4sW1AsJoVZQELkWC+usVvbt38/qX35h8/LlBAUH02/YMG5/6imPBgyr0chlXr4zP2mzlesf6eiBA/zw7rssnjkT2zmjRhq2aMGlV1/NoLFjadSyZdH2EIOBHhdgsFienl5in44oXd2gIGLDwgjx0HvmcH4+6zIzK/36flFRhFWyaUwptVZrnVDscxIo/qc6BYp1v//Ox08/zaFdu+jYowfDH36Yjj160K9/fyK3bj1r39S4OP789VcADmzfzs+ffMKir74iLDqae15+mZ7XXeexcgUbDHSLiPB4O26B08m2nBwO5OXhcDhY99tv7Nm8mVadOhHfvz9G1z9q6vHjfDJpEn/MnYvBYKDPzTcT37EjtRo1wlK/Pts3b2bNr7+y6c8/MRqN3P3KK1wxfHhRh/aFFixS7Xb+rCbv+UAVpBQdrVaaVLF2cTg/n/WZmVUaniyBogJqcqCw5eXx2eTJ/N9HH9GoZUueufVWblmwgJVz5+IMDaX2H39gysnBHhmJPSKC0JQUlMPBsWuuwZiVRbcRI9g3bhwrWrbkrQceIGnjRnoMGcJ9b76JJTzcI2U0KkVCeDj1PDQaKtVuZ31WFtkOBw6Hg0m33MKudevIz8kh2GLhovh4Js+Zw7ZVq5j297/TPS2NK/v0Ieq116jdsCH9unUjbO9eANI7diTluuvYdtllPPvyy2xZsYLeN97Iva++WvT7W4xGekZEeOwuMZCtycjw6ZKaNVkDs5kuYWGV6rvwRJAA7wUK6aOoRg7t2sXUcePYv3UrN48dy4tGI61feIHcJk0IPnGC3JgYTl1++Vmvybz44qKfQ44eJSgjg4Q776TFZZdx0SefMHPuXL566SUODx7MU199Rf2YmCqX88waBu0sFlqFhlZ6+KnWmt25uexyW6hl3W+/sWvdOvJcQznzsrPZuXYtb0ycSPB337HUbKaT3U7Orl381qABALsffBBzaiqmrCzqLllC+xdfJHT0aJ77/nu+nT6dmVOncuzgQZ7//nuCQ0PJcThYlZFBj8jIGp2QLsvH6y7XdEdtNlLT0uhotdIoOLhcr9Fasz8vj60emOho8OI6FVKjcJOWnc3WJUs4fdllXihV1axasIDX7rkHc2gok556ir+/9x7hu3ax/4472PbMMzjCwsp3IIeDmK+/psPTT+MIDWXtjBn85nTyytixmMxmnvziC9pdeqnHyl3PbCauEndZeQ4H67KyzhvN9M20acx85ZWz0je0AGYAVwJZTZvyy5Ah/J/FQkx8/FnNUmeEHD6McjjIjYkhfNs2Tn32Gbd8+indBg/mXx9/XLR/lMlEj8hIjDV0nsXGrCxJHe8ltYKCaG+xUKuUSZ0ZBQVsysryyLDk8O3b6Xrrraj33iN0yJBKHaO0GkXNvV2qhJBnn6X7jTdSa8UKfxeliNaa7956i5dGj6Zp27a8sWQJty5YQPDx46ycM4fN//53+YMEgNHIwb/9jT8WLcJWqxat33yT2N69mfrLL4SGhfHkDTeQuGiRx8p/3Gbjv2lpJOfllSs3j9aaA3l5LE1LK3bIa6tOnQg+ZzZ4PSA+OJiNzz9Pn5gY/vbFF3zx2mtMHTeOSbfcguOcET15jRuT66o5NfviC2769FNWderEmp9/5tNJk4r2S3P9I9dE+U6nX1ZKu1CctttZnp7OKtfa4KftdgqcTjIKCjiYl8eGzEyWpaV5JEhErV1Lj+uuQxUUoJuWmDe1SqRG4Sb9+HGMPXtiPnWKP375hRy3ETH+YLfZePfhh/lt5kx6DR3K/W++SbDFgvnECcxpaWS1aVOl4xuzsjDYbNhr1cKYnU1qbi6Thg/nwPbtPPrxx1w2cKCHfpNCESYTF1ss1C2h7+K03c6W7GzSS/nnOdNHsTsxkd65uSykcBTT+7//TuKKFUwdN66oWQogxGrlkRkzuPTqq4s/oNNJ6zfeoP2LL7K8RQv67tvHXa+8wqCxY4t26RIWRowXhkH6086cHFkLvSbQml4DB2I+dYqV337LZV26eKWPQmoUbnR0NH99/TUoRbdbbyUoLc1vZcnJyuK5227jt5kzGf7QQ0y79lp6jRuHIS8PW926VQ4SAI6wMOy1amHIy+OyYcPo9cILPD9nDi07deKVsWNZPm+eB36T/8koKGBVRga/nj5NYkYGe3Jz2Zuby7rMTH5LTWV5enqpQQLAaDQy6csvmR0ayi/AfUOG8O7KlaiwMPZs3kz+OR9++Tk57N28ueQDGgwkPfggW597jp779vFr3bp8+cwzpOzZU7TL5jKCV3XjcLWLi2pM68I+CaVY89lnLP/pJ3KbNfPa6SRQnCOnRQvWfPYZoYcOEXvffX4pw+mjR3niuuvY9Mcf3PfGGzzUsSMJd9+NOTUV5YUPLGdwMKd69KDZV1/R64kneG7WLNrExTF13DhW/PSTx8+X73RyxGZjW3Y2W7OzOZyfT045J3zpvDxq9+rFkNOn+f6aa7jy00+L+hSKa5YKtlho2alTmcfde889bHnhBdo3a4bFbOatBx/E6Uql7tSaxMxMCmpIavXk/HxsNeR3uRAZcnOJmzCB+LvvBq3Jb9CAfNfADa+d06tHr6ZOd+/O+nffZcfjvk8zdWjXLh4dPJiUPXt4+quvGBUVxSV//ztp8fGsmjWrYv0R5aUUO558ku1PPEGT776j2+TJTJ41i4vi45k6bhxrXHMw/E3l5VG/Z0/6JCfzZf/+BH311VnPx/fvz0Xx8YRYrSilCLFaucjVoV0e+8aNI/Gnnxjx3HNsXbmSDx57jG+mTWPNr7+SabOxMwDSPVeV1pq9NeD3uFBFbtzI5QMG0Pjbb0nv1Al81HUgw2NLcGSoa1VWrbEcOEBO8+ZeP+eWFSt4cfRoTEFBvPDjj/Q4dYpL7rqL9C5d+MtbQcJN0oMPYszN5aLXXye3aVMmffMNT990Ey/dcQdP/+c/xPXt69Xzl2XnxIkMOniQ9y+/nEYzZ573vNFoZPKcOWxfsoRDW7fSqUsXrhwwAIdShesSOJ2kFxSUmrZZm0wMuOIKrjOZ+Nenn7IAigLOlDlziAkOJtzPCRCr4rjd7rFFcJxOJ0f37+fgjh0Eh4ZSv1kz6jZpQlAAZROuKZTNxkXTptH6jTew1anD6q+/5viVV/rs/NX3He8jLd99l7b//jd//vILme3be+08y+bOZfrEiTRo1oxnZs6kQbNm2DZu5GTPnqz96CMKPDQZriw7H38cbTSSMnQo1ogIJs+ezZM33MALt9/OszNn0rFnT5+U41w/zZjBhz/+yNZBg7j+s8+KnZthNRppHx7O9SNGlHgcp9akFhRwwmbjYH5+sSv1bVyxgtscDuYA8UB6dja71q1j7W+/UW/wYLpHRnrwN/OtqtYmTJmZqFdeoeHXX9M4M5M+wDbgTtfXdmBzTAwZEybQ8dZbMdewQQD+EpSZSbMvv+TwTTex9fnnsUdH+/T80vRUhsM33URBWBiX3n47QadPe/z4TqeTmVOn8uq4cVwUH88rP/9clGgsvUsX/pozhwJffjApxa5//YvsVq1AaxodOcKUOXOo37QpU269lW1//eW7srjocePY8eSTXDZoENd9/HGxQeIii4W+UVE0LGOik0EpagcF0c5qpX90NBdbrefN8diyfz8jtKYp8IFr25lO8ZN2O0eq6bDSVLudk5XM6RSUmkqzsWPp06YNAz74AGNuLn+1a8ctTz7JqwsX0v+55zjctSst69blwYMHeeJf/yKoZUvmvPAC+dLUVXFOJ/UWLybunnvA4cBWuzZL//tfNrzzjs+DBMjw2LOUlMIjKjGRHkOHcrpbN/6aNQvtoZXRcrKymD5hAqvmz+eK4cO599VXqZOURPcbb2T3P//J3nvu8ch5KqvV22/T9qWXWPXttyS1aMGT11/P6WPHmDxnDu0Sih1F53EFjz3GDR9/zLcNGqDWrCn2DrVNaCjtrNbKn8PpZHVmZtG8jTW//srUceP4Z3Y2LwADgSUWC49++CGXXn01oQYD/aKjq91EvKqk61j55pv847nnWGQ2c2D0aNo/9dR5AwfOsOzaRei//41z5UouP36cejExTHzmGWLPNOdWU1prknfvZmdiIjvXruXgjh1kpaWRmZqKLS8Pa2QkYVFRRNWpQ0z79rTs1InWXbrQpE0bnE5nifnJ3IUeOEDT2bNpOnMmlkOHyG3cmBXff1/uofqS66kCvJHrqcmsWcRNnMi+sWPZ8sorVS0iB3fs4JW77uJwUhJjJ0/m2nHjiNy6le433ojDYmH5jz96dbhbeQSlptJr0CDMp07x5/z5HLRaefy668g4dYpnZ8/2erA4NXUqo/79b1ZGRHAkMZHQYu6kYkJC6OKBvpszaUdO2u1FczX2rV3LXzk5HAdGt2/P9CVLiv65O1ittAwNrfJ5fSWzoIClFRzurQoKaPrVV7yUksKs11+nV79+3Pvhh4RVoIa7eflyfvjnP/ll715+aN2a0LlzCfPyCB1PS05KYumcOSydM4eYQ4eoAzQOCaFFw4aYw8I4Xa8eW1u1Ijs9ndZ79nAyLY3tBw5w3G4nDYioXx+0Jis9HXt+PiFWK23j4nj5rbcIP3SI/Hr1yG7Viqi1a7l8wAC0Upzs04eDt97K0UGDcJYzHQhIoKgQbyUFvOiVV0jv3JljVZiI5nQ6mffBB3z5wgtYwsN5eMYMulx+OeFbt9LjxhtxhISw4ocfyGnRotLn8CTL/v30GjiQgrAw/lywgMP5+Tx5/fWkHT/O459/7rUO7n0ffsiYJ54gJTSUzStWENKkyXn7NDCbSQgP99hSpg6tWZORwQlXsFj3229k//473/78M47wcN5etgyj658wxGCgf3Q0hmpSq1ifmVmxmdgOB7Hjx9N07lyuAgyjRjF+6tSi378ijCkpRI0aRY/Nm9kUFMRfr71Gg1L6kQKB1ppDX39N3vTpHNu/n2kGA1169+b3deuIysg4a9/DQ4ey7qOPABjQogVB58zmn12vHsOPH8cAHAVCgFD+10Gc9I9/sP2ZZ1A2G80/+4yjAweSW8kZ1jUyUCilBgBvULhm9kda65fPeT4Y+AK4BDgFDNda7y/ruL7IHhty+DB5jRtX6PiH9+zh3YcfZvOff9J1wAAmvvYaUXXrYszO5oquXdEmU0AFiTOiEhPpccMNpMXGsmLePFKPH2fS8OEk79rFQ++/79E05QArfvqJRnfdxTCjkZW//46xXbvz9rEYjfSJjMTk4aR9Dq1ZlpZ21siglT//zLQxY3hoyhS6jx9ftL1zWJhXFq7xtByHg99TU8ufdM7hIO4f/6DJnDk8Bhx45BFGPvJIlQOy89136Tl5MmFOJ3OuuIKIr7/GEGAZekMOH8Y8ZQoxP/9MC1dg3dS6NevmzqV2gwZEr16NNhiw1apFgasmq81m7FFRAESvXk1QWhpB6emYU1MJSk/n56QkHvr+ewDeAvKBPOCwUhi6dKHWDTfQdvhwImrXrnL5a1ygUEoZgV3AVUAysAYYqbXe5rbPvUBnrfU9SqkRwA1a6+FlHdvbgaLu77/TddQo1r/zDik33FDm/qnHj/PNq6+y8IsvCA4J4c7nn+eq22476x+v8bffkhof7/e0ISWpv2ABzpAQTvTrB0BWejrP3XorOxMTGT1pEtePH1/lDxKtNXOmT+erF1+kfUICL732GoZiRpopoEdkZKkJ16rilN3OCrf3gXY6iWnZkrSCAk7v2UOQqynAYjTSLyoq4GsVm7Oyyj8T2+kk9r77aDprFk8Cm+66i3EvveSxsjj276fu0KFsTUlhRv/+/PO99wj3Q+dscQ7u3EnDYcO4PiWFP0JCODpwIBGPPUZBFf8nz/R5uaeWMYeEcMmVV7J382aOHTiAwWCg3aWXcuk11xDbpw/NO3Qotg+jLH4LFEqpRGAjsNn1tUlrfaJSJTn7uN2BZ7XW17gePw6gtX7JbZ+Frn1WKqVMFNbc6uoyCl3ZQHHLyJHYo6Np2akTLTt2pEmbNsVWtY3Z2XQbMYLoNWtY98EH/5tz4UZrzd5Nm1g8cya/f/MNtvx8rrn9doY/9BDR9eqB1rT84AOymzWrUlOWP0Rs3kxGx47k5+by+oQJrPi//+OywYO57403KtR+7c6Wl8cH993H8LlzWThoEMM/+KDEoZUXWSy0LaEj1VM2ZWVxwO3D1fTQQwz84gveHT+eplOmFG2PCwur8oI13pTncPBbWhrOct4QRm7cSPdBg3jJ4eDb7t2ZPHs2Jg8HZO10suizz3j/qafoVbs24//1L0JHjfLoOcpXEE3t5ctp/dJLvNyoEW/89BMXhYZy3bhxdHvoIY/NByltDRWDwcCeTZtY/csv/PXLL+zbsgWA0LAw2iUk0KJTJxq2aEGjli2JqluXEIuFYIsF7XSSnZ5OdmYmp48e5ci+faTs3Us9pfj8k08qVc6qBopGQGfXVwIwGDipta5ST6tS6mZggNb6LtfjvwHdtNYT3fbZ4ton2fV4j2ufk8UcbxwwDiAmJuaSAwcOVKg8drudHr16sWnTJmyuDwhzSAjNL76Ylp0707x9exq1akXj1q2p3bAhQTk5XDZsGLXWrCFpwgS2P/44p1NT2bVuHbvWrSNx8WL2b91KUHAwPa69lhEPP1y0VnNISgqdHn2UBgsXcmjYMDa8807FL6Cf1Fqxgh7XX8/Of/2L3Q89hNaaeR98wGeTJ1O3SRPumz69wnMttqxYwYePPsrrO3dyLbD2k084eu21xe4bbTLRMzLSY/0SJbE7nSxNSyPPNc9C5eUR17IlJw0G9u3fX3QDEWY00jcqyuvlqawtWVnsq0Bep9Tjx3mtTx+OhIby2qJFHmkOKcnOtWvpPHQo/fPz+eWqq3B+/jl4qZZ4rqi1a2k/ZQp1VqzgiMHAOKeTgtGjue2xx4isU8fj5zvT57V382ZaljLq6WRKCltXrWL7qlVsW72a5N27KSjnSDVrZCRdOnXiz2XLKvV+9GjTk1KqPXCz1vq5Cpfk7ON4NFC4q0rT09JTpziclMSezZvZu3kzezdtYu/mzWS7dWAZDAYsERHUDg9ncno6ozMyuNlk4jtXHiajyUSbuDj63nILvW+4gTBX+yUOB82++or2kydjKChgx2OPFQ6BrU6L47g1Tex4/HF2//OfAGxfvZpX776bE8nJJFx5JX976iladOhQ6qFOHTnCF88/z5+zZ/ONxcKNOTlseuUVDrhlbnVnUIq+UVE+W6b0SH4+iW7rF+c8+ijDP/2Ur8aPJ9KtVpEQHl7m/A1/qEhtIvTAAaJXr2bizz+zdvFiXl24sMy/nydkHziAZcgQrjt6lE21a3Pghx+gmD4pT+py//3EfP01p4KCeNZu5/dOnbhz2jTaxMV59byV4XA4OHn4MCl795KZmkpedjb5OTkYjEYsERFYIyKIrFOHRi1bEh4d7demp2Za6wPnbPtSa/23SpXmf8cIuKankvootNacPnaMw0lJHE5K4vTRo2SlpZGVno7T4SDOZuNwixZEN2jAiEOHCOvYkdODBqENBpTDgQ4KwmGxEL1mDb0GDeLE5Zezadq0gOu0LjeHozBYzJ7N9qeeIun++wHIz83l/z76iG/feIOcjAw6dO9O59696XL55dRu2BBbXh75OTlsX72a5T/9xLZVq4gwGlnauDGxBw6w7emn2VNKIkZ/DEldkZ5eNL/CYbPRrnlz9lks2HfvLrprqxUURM8AnK1d3tqEIS+PnoMHY05KIiYnhyFPP81NPkyI6XQ6Sf773xk1bx4GpfjujTeoPXKkR89hPnUKW61a5OfmkjtqFPuWL+ddq5WhTzzBwDvuqFR/QCDyZ6BYAcQA+yjso0gDBmqtqxR+XR/8u4D+wGEKO7Nv1VpvddtnAtDJrTP7Rq31sLKO7bc1s7WmT+/eROzYcdbmfXfcwZZ//xu0pt7ixYU5WgK0qaLcHA7iJkygyXffsWzxYtK7dCl6KistjXkzZrBm4UL2bt5c7IJFMe3a0ePaa7l64EBuePBB9t9xB4duu63E0/mqyelcZxagOeOvt9/mxcmTmTRrFvFXXFG0/fLISKJ81GxSHhWpTVw8aRKt3n2XEeHhrG3Zkqm//FKpYbBVdfD77wl+6CHGZ2Ux+J57GP3AA5iq2PQVdPo0LT78kFbvvssXt97KY7/8wonkZPoNG8aYSZMK+wtrEL+OelKF/52tgE5ALWDhmeagqlBKDQKmUzg89hOt9QtKqSlAotZ6nlIqBPgSiANOAyO01nvLOq7fAgWA00ntFSuIWrcObTSijUbyGjfmSAlt7tWZKiig4U8/FY38UjYb+pwOwIzTp9myfDnZGRmYQ0MxBwfTpE0bLjlwgNRLL8UeFVW4Mlcpb26DUvSOjPRbMr5V6emccNUq7DYbd196Ka0aN+bpuXOLJkM1Dg4m3kf5uMqjvCOdotesoefgwcxr0YKbDx3itcWLae62zrqv5WRl8fmUKaz89FO2G40ciYsj9d//JrscqeLdhaSk0PK992j2xReYcnJYEBHBfRkZ0KULY6dMoWOPHl76Dfyrxg2P9Sa/BooLVMSWLXS97TZ2338/ycOGlZjpNvjoUTo89RSNf/yRXQ89xM7HHivz2O0sFtp4eZRTac59X/zx4os88frrrLn/fvKeegooHLJ7ZXQ0IQHQhFHe2oSy2+nTpw+O1FSanjzJtY8+yshHHvFRKUu365dfsDzwALedOkUEsL1TJ3JGj+bE0KFFcxbOZczJweF6n1zeqxfhu3fzg9XKpMxMTjVvzt+eeIKeQ4diqE59ghUkgaICJFD4nnXPHmInTqRWYiL2sDCShw/ndNeuhbUNpYj58kuazJ5NrdWrcZpM7H7oIZImTjyvBnKuCJOJyyMj/T5XYXVGBsdco09yMjJo0aYNTUNCWJuUVJT7q6o5pzylIvMm6v3wA9MmTWJFaChvLlsWUCnCnU4n677+mvDnnmPE6dM0Ae66+mrqjhxJv+Rk2q9aBWFhhCYnY9m9G2N2Ng8+8giJK1ZQ67//ZXt+PuYuXbh23Dh6XX+9X3+3UIOBUKORfKeTHIej/JMfK8hbgULSjAuPyG7ViuXz5xO1bh0tPv6YmC+/JOY//yHlxhsBqPXXX5gyM9n94IMcGjas3BMLO1utfg8SAG0tlqJAYYmIYOGAATw9fz6RM2aQNmECAPvz8mhjsfg1WWCew8HB8qTqcDrBYODj48eZk5LC0//5T0AFCSgcXZgwahSOESP48PffSfnqK2atWEHWr7+STmHHabhS7NKancBO4OtnnyW0eXNqjxzJxGHDaJuQ4JehyyalaBIcTFPX+iXu7wmtNdkOB0dtNpLz88n00Pog3iQ1CjdSo/AcU0YGwSdOFKYrp/j+i7I0Dwmhk5cXa6oI9+yrxw8eJO6SS2gYGcnaHTuK+lg6Wa0092OywHKNdHI46HH99ey9+mquePNNWsfGMnn27ICdC+LObrORtGEDxw4e5PjBg5w+ehSLa4hoZJ06tEtIoIEPFhkrSbDBwEWhoTQJDi53epn0ggKScnNJ8UD6eqlRiGqlICKCgoiIoscVDRIhBgPt/NgvUZw2oaFFgaJeTAzfJiTwYmIiYYsXkzlgAABJubnEhIT4pRaU53BwoBwfNs2++ILaq1bxcWgoORkZ3DllSrUIEgBBZjPtu3alfdeu/i7KWQxK0SokhNahoRXOPxZpMnFJeDhtQ0PZlZvL4QBc76Tm9uqIaq2j1UpQgHU6RgUFUddtCGz05MnEAzMPHSralut0+u0ffU9eXpkd2OaTJ2n3wgskx8fz+H//yzWjR9PMiys3XgjqBAXRLyqKdlZrlZJUhplMxIeH0zsqitoBNNQaJFCIANTAbA7Imc7AWaOv2nXtSu6llzJvxgwcbmkWknJzi5074k35Tme5OrDbP/ccpuxsHgkLI9hi4dZHH/VB6Womk1J0Dguje2QkFg+Odos0megRGcml4eEePW5VSKAQAcWoFB0DYORQSWoHBZ2VtXboPfdwz/79xPXpA67gkOVwcKSSK8lV1r7c3DJrE6EHD9Jk1iwSb7iBb5Yt44YJE7yS1+hCEG0y0Scqyqtp5hsEB9M3KoqLLBa/D+iQQCECSjuLhdAAuYsqyUVundWXDRpERnQ0zZKSqLNsWdH23T5cJ9rudJYrVUduTAx/LlzI/UeOEFG7Ntf5eand6qp1aCg9PFyLKIlRKdq61oOv78dRaRIoRMCINJloEcApu8+oazYT5RpZYjSZyJw4kcNAk+efL9ono6CgaDitt+3Py6OgjNqEyZXUcllGBqv+/JNbHngASwCNKKsOgpSiW0QE7f0wZNtqNNI1IoJuERE+S4rpTgKFCAgK6BIWVm1G31zk1lfRb/Ro3gwKoumGDUS5Dcvemp1d7nUgKsuhdZm1CVNmJv26d6fVm2/y5QsvUKdRIwaOGePVctU0FqORXpGR1PPzXJN6ZjN9o6LoaLUS7MPBHhIoREBoERpKpJ9yOVVGfbO5qLxhkZHsGTmSU0CzqVOL9sl2OLzeBHUoL49815oZJWnz2msEnzjBEpOJXWvXMuLhh0tcFEqcr1ZQEJdHRlZ6foKnGZSiRWgoV0RF0c5iIcgHN1cSKC4AUSYT7S0WLg0Pp6PVSqvQ0KKmk0AQYjDQ1o+T1CrLva/iqnvv5XZgqmuC4RlJublkudYp8TSn1iSVEYise/bQ8oMPODhiBFNnzaJRy5b093AK75qsgdlM94gIzAE2VBvAZDDQxmLhqlq16BIWRrgXm6QC77cXHtMoOJj+0dFcHhVFa4uFBsHBtAgN5WKrlcujougeEUGdABiv3amK48/9pUFwMBGugNuoVStOXnMNX86dS77bh7dTaza7rZXsSYfz88ktozbR4emncQQH88Ull3Bg2zZGPvqoX1KIV0eNg4O5JDzc7yOOymJUipiQEPpGR3ut/6L6/XeKMpmUIi4sjEvKGIddx2yme2Qk3SMi/DZeu4HZTIMAnTNRHu61iqHjx1Pr5EnaX301YTt3Fm0/abdzqALLkZaH1rrMZq3QgwepvXw5O//5T2Z88AFN27al1/XXe7QcNVXT4GDiwsICPkicy1t9fBIoaphwo5E+UVE0qUAbdB1XB1mr0FB8+W9hCvA5E+XRwGwuqlV07NGDuh070mLnTlq//vpZ+23OzibDg01QKTYb2WUkk8uNiWHJqlV8Xrs2ybt3c+ujj9aYldy8qWlwcLUaWOELEihqkCjXKnCVqR0YleJiq5VekZFebet0d7HVGvBzJsqilCrqX1FK0ff++3lXaxp//z3WPXuK9nNozZrMTOxlNBWVh9aa3Tk5pe5j3bMHtCarTh2+ev11WnToQPchQ6p8bih8r1iNxmp3t10eDc1mCRLFkEBRQ9QOCqJ7RESV8yNFBQXROyqKNl6uXdQJCvLqrFZfahAcXDQCqvuQIXzWuDE2oPUbb5y1X47DwdrMzCqn9zhms5Wamjr46FF69+9P25df5vdZszi6fz+3Pf54lRbsOdOpe02tWgyqXZsroqMZUKsWPSIjuchi8elQTW+pZzYTHx4uQaIY1f+vK6gTFES3iAiPdQgblKKdq8PbG7ULk1J0qWGTvc5kujWaTPScMIEPtKbx7NlY9u07a78TdjvbyqgNlGVXGX0TF0+ZgrLb2Xv99XwzdSoXxcdz6dVXV+pczUJC6B8dzaUREdQxm88a/WNUitpBQbS1WOgfHc3FVqtPhmp6Q62gIBKqQce1v/glUCilaimlFimldru+R5ewn0MptcH1Nc/X5awOIkwmEsLDvbJYTqTJRG9XrhlPHv1iqzVgkp15Sj2zmWhXreLKW2/l7chIPmnZkvy6dc/bd29uLnsqOb/icH4+6aX0ddRZtowmc+awZ+JEvv3vfzmZksLfnnqqwnfJYUYjPSMj6RwWVq6/lVEpWoWGcmV0dLWrKUaYTHT10v9QTeGvGsVjwG9a6zbAb67HxcnVWse6vq7zXfGqhxCDgW7h4V5Nx21w5Zq5PCqKMA98uDc0m4mpxqOcSnOmVhFitRJ/553cnZTE/pSUYvfdlp1d4ZFQdqeTraUMtTXk5dHpkUfIbt6cTXfdxZzp04nt25cul19eofO0CQ2lT1TUWckPy8tkMBRmVI2IILQaNEdZjEav/w/VBP66OkOBz10/fw5c76dyVFsmpbgsIoIQH92Zn6ldVCUXU5TJRFwNbgOuYzYXrSNw3bhxBIeGsv/hh+n80EPF7r8xK4ujFVi7Ymt2dqmzsC0HDmDMy2PT1Kl89+mnZJw6xd+efLLcxzcbDFwWEUE7D+QyOjOSrmkA3xQEu35fX/0PVWf+ChT1tdZHXD8fBeqXsF+IUipRKbVKKXV9aQdUSo1z7Zt44sQJT5Y1IMWFhRHu44lTRqXoGBZW+M9VwTswiyupWU2v3rd11Soiatdm8F13kb9yJc2++IK6S5act68GEjMzSS5HzeKkzcahMoJKVtu2/P7XX+zp2JEf332XHtdeS5vY2HKVu1ZQEH0iI6nrwVxGJoOB2PBw4sPDMQXY391sMNDdTwn2qiOvBQql1GKl1JZivoa676cLh4CUNAykmWsN11uB6UqpViXsh9Z6htY6QWudULeYduGapGVoqF8nqdU1m+kTFUXDcn6omA0GuoaH14iRMWWpHRRUVKu4YcIEPrdYSLFYaD9lChRTG9DA+qysUvssHFqzqbTZ3Q4HzT7/HGWz4QwJYfbrr5Ofm8uoxx8vV5mbBAfT3Yt31o2Dg+kTFRUwaWOCXLVxX99oVWde+8/VWl+pte5YzNePwDGlVEMA1/fjJRzjsOv7XmApEOet8lYXZ/I2+ZvZYCAhIoL48PBS78oaBwfTLyrqgvqnLKpV1KrF1ePG8XBODpFbttDs889LfM227Gw2Z2WRd86w15M2G3+kpZU6ua7VO+/Q+eGHqb9oEYd27WL+J59w1ahRNGnTplxljfPBaB+Lq3O8uZ87uk2uVOHVKQFlIPDXLd48YLTr59HAj+fuoJSKVkoFu36uA/QEtvmshAEoSKmAyz1zJhAkhIdTKygIq9FIpMlUNGQ3Pjw8IBOqeZN7reL6e+/lB6uV1XXqcPHkyQQfPVri6/bn5bE4NZU1GRmk5OeTmJHByoyMUudMRGzaRLuXXyZlyBCODBzIR08+SYjFUq7aRJewsLPSpXubQSk6hYUR76cRRkFK0T0igugAyG9W3fgrrL4MzFZK3QkcAIYBKKUSgHu01ncB7YEPlFJOCgPay1rrCzpQdAzQYaVKKRoGBwfsOtf+0NZiYUV6OuHR0Vx7990Me+01Pnj4YWz1S+qOK6SBozYbR8ux6JEhN5f4e+/FVqsWm6ZNY82iRaxfupS7nn++zCVOO1qtxPjp7r5xcDARRiNrMjPLTEPiKSGujusLqWbrScrXi8D7QkJCgk50W0CmvNLsdv5IT/dCiaqublAQl0VG+rsYogJWpKdzym4nKz2d8d2707B5c17+v//D6HCgPXBXe/GkSbR6911WzZpFSs+eTOjViyCzmTeWLsVUyvHbWiw+rUmUxO50si4ri+NeXgnQYjRymXRcl0kptdbVJ3yeC6tNoJoy1sCZzBeCNq4cUGGRkYx5+ml2rFnDiaee4oquXQk5fLjKxz80fDjbn3ySE1dcwY/vv8/R/fu587nnSg0SLUJCAiJIAAS5Bjm08eJaJHVciw5JkKgaCRTVQHuLpdonz7sQua+t3W/4cNomJPDet98SlJZGt5EjCUpLq9RxQw8cAK3JvPhikh54gIM7dzJz6lQuGzSI+CuuKPF1tYOC6BBg2XqVK12MN7ILtAgJoVuALjpU3cgVDHDRJpPfR4qIyjtzt2wwGLj75ZdZn5bG5L59CUtK4tLbb8dQwdnZEZs20bd3b1q+9x4ABXY7r997L6FhYYx3W4b1XGaDgfgAzoraMDjYY3f+Qa71WDpWw/UkApUEigCmgM4B/M8tylbfbC5KrNi6Sxeuuf12XlqwgAWPPkrtlSuJGz8eytmhW2fpUrrfcgu2WrU4fOONAMyaNo09mzYxYdo0ouvVK/G18WFhAT8DOdxkok9UVJVujOoGBdG3guuxiLJJoAhgzUJCihbFEdWTUoo2bn0Co554glr16zP2yy9Z+/jj5DVsCOfcCDgcDtb8+ivfTJvGml9/xWG30+a117hs2DDy6tdn5Xffkd+gATvXrmXO9OlcMXw43QcPLrEMbUJDPTrj2puMriG03So4+99iNBIbFsZlkZEBHxCrI/kUClBmg6EoyZyo3hqZzew0Gsl2OIioVYtHP/6YJ667jnvWrOGpr77CYDAQvXo1tZcv58CwYTw1cSK71q8nPzubYKuVQW3bMmfTJg7feCObpk3DYbVy6sgRpv7979Ru2JC/v/hiieeONpmKJgBWJ/XMZvpHR3M4P5+9eXklrg4YZjTSOjSUxsHB0szkRRIoAlR7i0UyWtYQSikuCg1lfVYWAO0SErjzuef44LHHmDN9OsMfeogGv/xC67feov2LL3IZsBc4DVyVnc38nTt577nniLnzTlCKzNRUJg0bRubp07zw449YIyKKPa9RqWqdhNGgFE1DQmgaEkKa3U6O00mu04nd6STKZCI6KOiCSAsTCCRQBKAokymgs26KimscHMzevLyitSQGjR3LjsREvn7lFeo0bgzPPMORQYPY89pr5CxaRGtgC4X9VPk5OaxMTydGKfKys5ly662k7N3Ls7Nm0bpLlxLP2dFqrTHDQqOCgojydyEuYBIoAowCWbO3BlJKcbHFwsqMjKLHE159ldNHj/LGP/5B8q5djHrySfaOGcPUFSvIc0sCGGKx0LJTJ04kJ/Pavfeye906Hv34Yzr36lXi+RqYzX6beS1qHqm3BZjWoaHSgV1D1TGbqefWqRxitTJ59mwGjBnDd2+9xUujRxPTvj0XxccTYrWilCLEaqVNXBzHDh5kYq9e7Nm4kQfeeYceQ4aUeJ5Qg0EmaAqPkhQebvydwsNqNNInKqrGr9lwIcssKOC/aWln5dXXWjP/k0/48MkncToctOvalaYXXcSpw4dxOp2cPnaMgzt2ENu3LxOmTaN+TEyJxzcqRa/ISLnZEBVWWgoPeTcFkC5hYRIkarhwk4mYkBAOuE20U0ox+M47ievXjz/mzuXPefNY9NVXAETWrUtM27bc/9ZbXDF8eJlNkvFhYRIkhMdJjcKNP2sUzUNC6CTNBRcEp9aszMjgtN1e4j7HDh4k1Golonbtch+3vcVC62o4FFYEBkkKGODqBmAOHuE9BqW4tIwFn+rHxJQ7SBiUopPVKkFCeI0ECj+LMJlICLDFiIT3mQ0GjySssxqNXB4ZSXMvZmAVQhoz/SjUYKBbeDgmmTR0QbIajfSKjGRvbi7J+fkUlLMZ2GwwUNu1imCT4GB5/wivk0DhJw3NZjpYrZKX5gJnNRrpFBZGe4uFozYbDgqzn575MhsMBCmFonD1O+16XubZCF+SQOFjVqORTlZrtUnSJnzDZDBIxlMRsPxSZ1VK3aKU2qqUcrrWyS5pvwFKqZ1KqSSl1GO+LKOnWIxGGgUH08lqpW9UFFdER0uQEEJUK/6qUWwBbgQ+KGkHpZQReAe4CkgG1iil5mmtt/mmiJUTajDQMDiYukFBRJlMsrqWEKLa80ug0FpvB8pqZ+0KJGmt97r2/QYYCgRcoIg0mWhgNtPAbJbJTkKIGieQP9UaA4fcHicD3UraWSk1DhgHEFNKioOKMhsMhBgMmF0dikZXB2OY0YjVYCDKZJIOaSFEjea1QKGUWgw0KOapJ7XWP3r6fFrrGcAMKJyZXZVjRZlMXGSxUCcoSFJqCCEueF4LFFrrK6t4iMNAU7fHTVzbvObMJKh60tkshBBFArmndQ3QRinVQillBkYA87x5QovRKEFCCCHO4a/hsTcopZKB7sDPSqmFru2NlFLzAbTWBcBEYCGwHZittd7qj/IKIcSFzF+jnuYCc4vZngIMcns8H5jvw6IJIYQ4RyA3PQkhhAgAEiiEEEKUSgKFEEKIUkmgEEIIUSoJFEIIIUolgUIIIUSpJFAIIYQolQQKIYQQpZJAIYQQolQSKIQQQpRKAoUQQohSSaAQQghRKgkUQgghSiWBQgghRKkkUAghhCiVBAohhBClkkAhhBCiVBIohBBClMpfa2bfopTaqpRyKqUSStlvv1Jqs1Jqg1Iq0ZdlFEIIUcgva2YDW4AbgQ/KsW8/rfVJL5dHCCFECfwSKLTW2wGUUv44vRBCiAoI9D4KDfyqlFqrlBrn78IIIcSFyGs1CqXUYqBBMU89qbX+sZyH6aW1PqyUqgcsUkrt0FovK+F844BxADExMZUqsxBCiPN5LVBora/0wDEOu74fV0rNBboCxQYKrfUMYAZAQkKCruq5hRBCFArYpiellFUpFX7mZ+BqCjvBhRBC+JC/hsfeoJRKBroDPyulFrq2N1JKzXftVh/4Uym1EVgN/Ky1/sUf5RVCiAuZv0Y9zQXmFrM9BRjk+nkv0MXHRRNCCHGOgG16EkIIERgkUAghhCiVBAohhBClkkAhhBCiVBIohBBClEoChRBCiFJJoBBCCFEqCRRCCCFKJYFCCCFEqSRQCCGEKJUECiGEEKWSQCGEEKJUEiiEEEKUSgKFEEKIUkmgEEIIUSoJFEIIIUolgUIIIUSpJFAIIYQolQQKIYQQpfJLoFBKTVVK7VBKbVJKzVVKRZWw3wCl1E6lVJJS6jEfF1MIIQT+q1EsAjpqrTsDu4DHz91BKWUE3gEGAhcDI5VSF/u0lEIIIfwTKLTWv2qtC1wPVwFNitmtK5Cktd6rtbYB3wBDfVVGIYQQhUz+LgAwFphVzPbGwCG3x8lAt5IOopQaB4xzPcxSSu2sZHnqACcr+VpvknJVjJSrYqRcFVMTy9WspCe8FiiUUouBBsU89aTW+kfXPk8CBcB/qno+rfUMYEZVj6OUStRaJ1T1OJ4m5aoYKVfFSLkq5kIrl9cChdb6ytKeV0qNAYYA/bXWuphdDgNN3R43cW0TQgjhQ/4a9TQAeBS4TmudU8Jua4A2SqkWSikzMAKY56syCiGEKOSvUU9vA+HAIqXUBqXU+wBKqUZKqfkArs7uicBCYDswW2u91Qdlq3LzlZdIuSpGylUxUq6KuaDKpYpv9RFCCCEKycxsIYQQpZJAIYQQolQXZKBQSt2ilNqqlHIqpUocSlZSChFXB/tfru2zXJ3tnihXLaXUIqXUbtf36GL26efq1znzlaeUut713GdKqX1uz8X6qlyu/Rxu557ntt2f1ytWKbXS9ffepJQa7vacR69XWSlnlFLBrt8/yXU9mrs997hr+06l1DVVKUclyvVPpdQ21/X5TSnVzO25Yv+mPirXGKXUCbfz3+X23GjX3323Umq0j8v1uluZdiml0tye88r1Ukp9opQ6rpTaUsLzSin1pqvMm5RS8W7PVf1aaa0vuC+gPdAWWAoklLCPEdgDtATMwEbgYtdzs4ERrp/fB8Z7qFz/Bh5z/fwY8EoZ+9cCTgMW1+PPgJu9cL3KVS4gq4TtfrtewEVAG9fPjYAjQJSnr1dp7xe3fe4F3nf9PAKY5fr5Ytf+wUAL13GMPixXP7f30Pgz5Srtb+qjco0B3i7mtbWAva7v0a6fo31VrnP2/wfwiQ+uV28gHthSwvODgAWAAi4D/vLktbogaxRa6+1a67JmbhebQkQppYArgG9d+30OXO+hog11Ha+8x70ZWKBLHmLsKRUtVxF/Xy+t9S6t9W7XzynAcaCuh87vrjwpZ9zL+y3Q33V9hgLfaK3ztdb7gCTX8XxSLq31Erf3UEkpdTytKil6rgEWaa1Pa61TKcwdN8BP5RoJzPTQuUuktV5G4U1hSYYCX+hCq4AopVRDPHStLshAUU7FpRBpDNQG0vT/clWd2e4J9bXWR1w/HwXql7H/CM5/k77gqnq+rpQK9nG5QpRSiUqpVWeawwig66WU6krhXeIet82eul4lvV+K3cd1PdIpvD7lea03y+XuTgrvTM8o7m/qy3Ld5Pr7fKuUOjMBNyCul6uJrgXwu9tmb12vspRUbo9cq0DI9eQVqhwpRPyhtHK5P9Baa6VUiWOXXXcLnSicZ3LG4xR+YJopHE/9L2CKD8vVTGt9WCnVEvhdKbWZwg/DSvPw9foSGK21dro2V/p61URKqVFAAtDHbfN5f1Ot9Z7ij+BxPwEztdb5Sqm7KayNXeGjc5fHCOBbrbXDbZs/r5fX1NhAoctIIVIOJaUQOUVhtc7kuiusUGqR0sqllDqmlGqotT7i+mA7XsqhhgFztdZ2t2OfubvOV0p9Cjzsy3JprQ+7vu9VSi0F4oDv8PP1UkpFAD9TeJOwyu3Ylb5exShPypkz+yQrpUxAJIXvJ2+mqynXsZVSV1IYfPtorfPPbC/hb+qJD74yy6W1PuX28CMK+6TOvLbvOa9d6oEylatcbkYAE9w3ePF6laWkcnvkWknTU8mKTSGiC3uIllDYPwAwGvBUDWWe63jlOe55baOuD8sz/QLXA8WOkPBGuZRS0WeabpRSdYCewDZ/Xy/X324uhe23357znCevV3lSzriX92bgd9f1mQeMUIWjoloAbYDVVShLhcqllIoDPqAwpc5xt+3F/k19WK6Gbg+vozBDAxTWoq92lS8auJqza9ZeLZerbO0o7Bxe6bbNm9erLPOA212jny4D0l03Qp65Vt7ooQ/0L+AGCtvq8oFjwELX9kbAfLf9BlG4sNIeCu9Gz2xvSeE/chIwBwj2ULlqA78Bu4HFQC3X9gTgI7f9mlN4p2A45/W/A5sp/MD7CgjzVbmAHq5zb3R9vzMQrhcwCrADG9y+Yr1xvYp7v1DYlHWd6+cQ1++f5LoeLd1e+6TrdTuBgR5+v5dVrsWu/4Mz12deWX9TH5XrJWCr6/xLgHZurx3ruo5JwB2+LJfr8bPAy+e8zmvXi8KbwiOu93IyhX1J9wD3uJ5XFC70tsd17gS311b5WkkKDyGEEKWSpichhBClkkAhhBCiVBIohBBClEoChRBCiFJJoBBCCFEqCRRCCCFKJYFCCCFEqSRQCOEDSqklSqmrXD8/r5R6y99lEqK8amyuJyECzCRgilKqHoX5f67zc3mEKDeZmS2Ejyil/guEAX211pn+Lo8Q5SVNT0L4gFKqE9AQsEmQENWNBAohvMyVBfU/FK5ClqWU8tRqbEL4hAQKIbxIKWUBvgce0lpvB56jsL9CiGpD+iiEEEKUSmoUQgghSiWBQgghRKkkUAghhCiVBAohhBClkkAhhBCiVBIohBBClEoChRBCiFL9P1dpibiJLygRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUJklEQVR4nO3dfZBdd33f8fenFrZrTLBsKYqwjSVPDKkJFDyq6wINYLuDwQGb4qGmAQRVR014CJSkWMadgWaaiU0fSNy0MCoGlAT8gAOxCSWN8cPkgVhUBsePGMtPWIosLQHRAI3B8O0f9yy5Xt+Vdu+9e3f35/drZueeh98557vnrj767e+cezZVhSSpLX9vsQuQJI2f4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMdy0JSe5M8tLFrmMhJXl2kluT/E2SXx7TPh9O8oJx7EttMdy14JI8mOTMGcvenOTPpuer6jlVddN897PMvAe4saqeVlWXjrqzJCuBtcDdI1em5hjuEpBkxQQOcwJw53w3OkBtzwV2VtXfjlSVmmS4a0no75UnuSDJ7m744p4kZyT5XeCZwGeTfCfJe7q2/yDJTUn2d0M7r+7b5ylJvtLt51NJrkzyH2cc84IktwHfTbIiyZYk93Xb3JXkNTPa/7sktyX5bpLLkqxJ8vmu/Re63vSg7+8G4GXAb3f1P+sgtT+htgG7fR5wR9f+iCSfTPLpJEcO/06oFYa7lpQkzwbeDvyjqnoa8HLgwap6I/B14FVVdWRVfSDJU4DPAn8M/CTwDuAT3dj2ocBngI8DRwOXA695wgHh9cDZwFFV9RhwH/BPgacD/wH4vSRr+9q/FvhnwLOAVwGfB94LrKb372ngWHpVnQ78KfD2qjoSeGC22g9Q20zPBW5Psh74c+Ae4LVV9Z1BNejJZRK/ikoAf5CkP6AOBb48oN0PgcOAk5NMVdWDB9jnacCRwMVV9SPghiR/SC8Ub6D3831p9R59+ukkXxqwj0ur6uHpmar6VN+6K5NcCJwKXNMt+29VtRcgyZ8C+6rqK938Z4AzDlDvXGt//6DaBngeUMCNwDur6poDtNWTjD13Tcq5VXXU9Bfw1kGNqmon8C56AbcvyRVJnjHLPp8BPNyF47SHgGO7dbvr8c+0HhSUj1uW5E3dHS37k+wHfhZY1ddkb9/0/xswP9chkQPVfqB6p+tMV9trgA/NFuxJ/Df+JOUbryWnqj5ZVS+mdwGygEumV81o+lfA8TMC7JnAbmAPcGwXgtOOH3S46YkkJwD/k96w0DHdf0J3ABmw3agOVPsTahtgffd6JvArSTb0r0xyS5IP0/t+9CRkuGtJ6cbLT09yGPC39HrD073bvcCJfc23A98D3pPkKd198q8CrgD+gt4Qz9u7C6Xn0BteOZCn0gvUqa6Wt9DrHS+EA9U+F88Dbquq24HNwGemrw0kWUVvHP+9VbVp3IVreTDctdQcBlwMfAN4hF5IXdit+w3g33dDJr9aVd+nF4iv6Nr/D+BNVfXVbt0/BzYB+4E3AH8IPDrbgavqLuC/0PuPYS+9C5Z/Pu5vsDvWrLXPcRfPBW7r9vUHwFZ61zUOpxf8n6yqb467bi0f8c/s6ckiyXbgw1X1scWuZSEleRewq6quXuxatHjsuatZSV6S5Ke6YZmN9Hq0f7TYdU3Ac4FbF7sILS5vhVTLng1cRW8s/X7gvKras7glLTzH2QUOy0hSkxyWkaQGLYlhmVWrVtW6desWuwxJWlZuueWWb1TV6kHrlkS4r1u3jh07dix2GZK0rCR5aLZ1DstIUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDlsQnVKUWrNvyuTm1e/Disxe4EsmeuyQ1yXCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lq0EHDPclHk+xLckffsqOTXJfk3u51Zbc8SS5NsjPJbUlOWcjiJUmDzaXn/nHgrBnLtgDXV9VJwPXdPMArgJO6r83Ah8ZTpiRpPg4a7lX1J8A3Zyw+B9jWTW8Dzu1b/jvVczNwVJK1Y6pVkjRHw465r6mqPd30I8CabvpY4OG+dru6ZZKkCRr5gmpVFVDz3S7J5iQ7kuyYmpoatQxJUp9hw33v9HBL97qvW74bOL6v3XHdsieoqq1VtaGqNqxevXrIMiRJgwwb7tcCG7vpjcA1fcvf1N01cxrw7b7hG0nShKw4WIMklwMvBVYl2QW8D7gYuCrJJuAh4HVd8/8FvBLYCXwPeMsC1CxJOoiDhntVvX6WVWcMaFvA20YtSpI0Gj+hKkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoJHCPcm/TXJnkjuSXJ7k8CTrk2xPsjPJlUkOHVexkqS5GTrckxwL/DKwoap+FjgEOB+4BPhgVf008C1g0zgKlSTN3ajDMiuAv59kBXAEsAc4Hbi6W78NOHfEY0iS5mnocK+q3cB/Br5OL9S/DdwC7K+qx7pmu4BjB22fZHOSHUl2TE1NDVuGJGmAUYZlVgLnAOuBZwBPBc6a6/ZVtbWqNlTVhtWrVw9bhiRpgFGGZc4EHqiqqar6AfBp4EXAUd0wDcBxwO4Ra5QkzdMo4f514LQkRyQJcAZwF3AjcF7XZiNwzWglSpLma5Qx9+30Lpx+Gbi929dW4ALg3Ul2AscAl42hTknSPKw4eJPZVdX7gPfNWHw/cOoo+5UkjcZPqEpSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNWikcE9yVJKrk3w1yd1J/kmSo5Ncl+Te7nXluIqVJM3NihG3/y3gj6rqvCSHAkcA7wWur6qLk2wBtgAXjHgcaVGs2/K5xS5BGsrQPfckTwd+DrgMoKq+X1X7gXOAbV2zbcC5o5UoSZqvUYZl1gNTwMeSfCXJR5I8FVhTVXu6No8AawZtnGRzkh1JdkxNTY1QhiRpplHCfQVwCvChqnoB8F16QzA/VlUF1KCNq2prVW2oqg2rV68eoQxJ0kyjhPsuYFdVbe/mr6YX9nuTrAXoXveNVqIkab6GDveqegR4OMmzu0VnAHcB1wIbu2UbgWtGqlCSNG+j3i3zDuAT3Z0y9wNvofcfxlVJNgEPAa8b8RiSpHkaKdyr6lZgw4BVZ4yyX6llc7298sGLz17gStQyP6EqSQ0y3CWpQYa7JDXIcJekBhnuktSgUW+FlKSx8U6i8bHnLkkNMtwlqUGGuyQ1yDF3aYly/FmjsOcuSQ0y3CWpQYa7JDXIcJekBhnuktQg75bRsuCdI9L82HOXpAYZ7pLUIIdl9KQ012Eeabmy5y5JDTLcJalBhrskNchwl6QGGe6S1CDvlpG07CzW3U7L6UNy9twlqUGGuyQ1yHCXpAaNPOae5BBgB7C7qn4+yXrgCuAY4BbgjVX1/VGPI2kwH6qmQcbRc38ncHff/CXAB6vqp4FvAZvGcAxJ0jyM1HNPchxwNvDrwLuTBDgd+Jddk23A+4EPjXIcSaOzhz9Zi32+R+25/ybwHuBH3fwxwP6qeqyb3wUcO2jDJJuT7EiyY2pqasQyJEn9hg73JD8P7KuqW4bZvqq2VtWGqtqwevXqYcuQJA0wyrDMi4BXJ3klcDjwE8BvAUclWdH13o8Ddo9epiRpPobuuVfVhVV1XFWtA84HbqiqXwBuBM7rmm0Erhm5SknSvCzEfe4X0Lu4upPeGPxlC3AMSdIBjOXZMlV1E3BTN30/cOo49itJGo6fUJWkBhnuktQgH/mrpviHr6Uee+6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3ybhnN6w4THwcrLQ/23CWpQfbcG+Y939KTlz13SWqQ4S5JDXJYRovKoSNpYdhzl6QGGe6S1CDDXZIaZLhLUoMMd0lqkHfLLEPeYSLpYOy5S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIO9z17zM9R57/xzf8uV73AbDXdKC84N3kzf0sEyS45PcmOSuJHcmeWe3/Ogk1yW5t3tdOb5yJUlzMcqY+2PAr1TVycBpwNuSnAxsAa6vqpOA67t5SdIEDR3uVbWnqr7cTf8NcDdwLHAOsK1rtg04d8QaJUnzNJa7ZZKsA14AbAfWVNWebtUjwJpZttmcZEeSHVNTU+MoQ5LUGTnckxwJ/D7wrqr6v/3rqqqAGrRdVW2tqg1VtWH16tWjliFJ6jNSuCd5Cr1g/0RVfbpbvDfJ2m79WmDfaCVKkuZrlLtlAlwG3F1V/7Vv1bXAxm56I3DN8OVJkoYxyn3uLwLeCNye5NZu2XuBi4GrkmwCHgJeN1KFkpYk711f2oYO96r6MyCzrD5j2P1Kkkbns2UkqUGGuyQ1yGfLTIAPYpI0afbcJalBhrskNchhmSWkpVvLWvpepOXInrskNchwl6QGGe6S1CDDXZIaZLhLUoO8W0aS5mg53QVmz12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhrkJ1RHsJw+rSbpyWXZh7t/n1SSnshhGUlq0LLvuS8Eh1skLXf23CWpQYa7JDXIcJekBhnuktSgBQn3JGcluSfJziRbFuIYkqTZjT3ckxwC/HfgFcDJwOuTnDzu40iSZrcQPfdTgZ1VdX9VfR+4AjhnAY4jSZrFQtznfizwcN/8LuAfz2yUZDOwuZv9TpJ7FqCWvzveJawCvrGQxxiBtQ3H2oZjbcMbe325ZKTNT5htxaJ9iKmqtgJbJ3W8JDuqasOkjjcf1jYcaxuOtQ1vqdfXbyGGZXYDx/fNH9ctkyRNyEKE+/8BTkqyPsmhwPnAtQtwHEnSLMY+LFNVjyV5O/C/gUOAj1bVneM+zhAmNgQ0BGsbjrUNx9qGt9Tr+7FU1WLXIEkaMz+hKkkNMtwlqUFNhXuSo5Ncl+Te7nXlLO1+mOTW7uvavuXrk2zvHptwZXdBeGK1JXl+kr9IcmeS25L8i751H0/yQF/dzx9DTQd8TESSw7rzsLM7L+v61l3YLb8nyctHrWWI2t6d5K7uPF2f5IS+dQPf3wnW9uYkU301/Ou+dRu7n4F7k2xchNo+2FfX15Ls71u30Ofto0n2JbljlvVJcmlX+21JTulbt9Dn7WC1/UJX0+1JvpjkH/ate7BbfmuSHeOubWhV1cwX8AFgSze9BbhklnbfmWX5VcD53fSHgV+aZG3As4CTuulnAHuAo7r5jwPnjbGeQ4D7gBOBQ4G/BE6e0eatwIe76fOBK7vpk7v2hwHru/0cMuHaXgYc0U3/0nRtB3p/J1jbm4HfHrDt0cD93evKbnrlJGub0f4d9G54WPDz1u3/54BTgDtmWf9K4PNAgNOA7ZM4b3Os7YXTx6T3aJXtfeseBFYt5Lkb5qupnju9xxxs66a3AefOdcMkAU4Hrh5m+3HUVlVfq6p7u+m/AvYBq8dYQ7+5PCaiv+argTO683QOcEVVPVpVDwA7u/1NrLaqurGqvtfN3kzv8xSTMMrjNV4OXFdV36yqbwHXAWctYm2vBy4f4/EPqKr+BPjmAZqcA/xO9dwMHJVkLQt/3g5aW1V9sTs2TPbnbWithfuaqtrTTT8CrJml3eFJdiS5Ocm53bJjgP1V9Vg3v4veoxQmXRsASU6l1/u6r2/xr3e/Gn4wyWEj1jPoMREzv98ft+nOy7fpnae5bLvQtfXbRK/HN23Q+zvp2l7bvVdXJ5n+UN+SOW/dMNZ64Ia+xQt53uZitvoX+rzN18yftwL+OMkt6T1WZUlYdn9DNckXgJ8asOqi/pmqqiSz3ed5QlXtTnIicEOS2+kF11Koja638rvAxqr6Ubf4Qnr/KRxK717bC4BfG7Xm5S7JG4ANwEv6Fj/h/a2q+wbvYUF8Fri8qh5N8m/o/fZz+gSPPxfnA1dX1Q/7li32eVvykryMXri/uG/xi7vz9pPAdUm+2v0msKiWXbhX1ZmzrUuyN8naqtrTBeS+Wfaxu3u9P8lNwAuA36f3a+CKrpc678cmjKO2JD8BfA64qPvVdHrf073+R5N8DPjV+dQ2wFweEzHdZleSFcDTgb+e47YLXRtJzqT3H+dLqurR6eWzvL/jCqmD1lZVf903+xF611umt33pjG1vGlNdc6qtz/nA2/oXLPB5m4vZ6l/o8zYnSZ5H7/18Rf973Hfe9iX5DL3hsUUP90Uf9B/nF/CfePxFyw8MaLMSOKybXgXcS3fRCfgUj7+g+tYJ13YocD3wrgHr1navAX4TuHjEelbQuzC1nr+7+PacGW3exuMvqF7VTT+Hx19QvZ/xXlCdS23TwXPSXN/fCda2tm/6NcDN3fTRwANdjSu76aMnWVvX7mfoXQTMpM5b33HWMftFy7N5/AXVL03ivM2xtmfSu7b0whnLnwo8rW/6i8BZ465tqO9nsQsY85tzTBeO9wJfmP4BoPdr+0e66RcCt3c/+LcDm/q2PxH4Uvcmfmr6h32Ctb0B+AFwa9/X87t1N3T13gH8HnDkGGp6JfC1LiQv6pb9GvDqbvrw7jzs7M7LiX3bXtRtdw+9nsy438uD1fYFYG/febr2YO/vBGv7DeDOroYbgZ/p2/ZfdedzJ/CWSdfWzb+fGZ2DCZ23y+ndAfYDeuPmm4BfBH6xWx96f+jnvq6GDRM8bwer7SPAt/p+3nZ0y0/sztlfdu/5ReOubdgvHz8gSQ1q7W4ZSRKGuyQ1yXCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQ/weac0V0A1ILOwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.57273275 0.44071084\n"
     ]
    }
   ],
   "source": [
    "u_pred, k_pred = model.predict(x_test, samples, [process_u, process_k])\n",
    "plot1d(x_u_train, u_train, x_test, u_test, u_pred[..., 0], ylabel=\"$u$\", ylim=[-2.0, 2.0], title=\"$N(0, 1^2)$\")\n",
    "plt.hist(k_pred, bins=30)\n",
    "plt.title(\"Histogram for $k_r$\")\n",
    "plt.show()\n",
    "print(np.mean(k_pred), np.std(k_pred))\n",
    "\n",
    "# sio.savemat(\n",
    "#     \"Normal.mat\",\n",
    "#     {\n",
    "#         \"x_u_train\": x_u_train, \"u_train\": u_train,\n",
    "#         \"x_f_train\": x_f_train, \"f_train\": f_train,\n",
    "#         \"x_test\": x_test, \"u_test\": u_test, \"f_test\": f_test,\n",
    "#         \"u_pred\": u_pred, \"f_pred\": f_pred, \"k_pred\": k_pred,\n",
    "#         \"noise\": noise,\n",
    "#     }\n",
    "# )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "feebcd9c",
   "metadata": {},
   "source": [
    "#### Half Normal prior for $k_r$\n",
    "\n",
    "$X\\sim HalfN(0, 1^2)$ is equivalent to there exists $Z$ such that $X=|Z|, Z\\sim N(0, 1^2)$. Choosing this prior for $k_r$ indicates that we at least known $k_r\\geq 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "960b478f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow_probability as tfp\n",
    "\n",
    "class HalfNormal(neuq.variables._Samplable):\n",
    "    def __init__(self, sigma=1.0, shape=[], initializer=None):\n",
    "        super().__init__()\n",
    "\n",
    "        self._num_tensors = 1\n",
    "        # initial value should be set to a non-negative value\n",
    "        self._initial_values = [tf.constant(0.1, tf.float32)]\n",
    "        if shape == []:\n",
    "            # make sure constant has at least 1 dimension\n",
    "            self._initial_values = [self.initial_values[0][None, ...]]\n",
    "        self.dist = tfp.distributions.HalfNormal(scale=sigma)\n",
    "\n",
    "    def log_prob(self, samples):\n",
    "        # Note: here, because a constant is considered, `samples` is a list of only\n",
    "        # one element.\n",
    "        return self.dist.log_prob(samples[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "19e7ac8d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Supporting backend tensorflow.compat.v1\n",
      "\n",
      "Compiling a MCMC method\n",
      "\n",
      "sampling from posterior distribution ...\n",
      "\n",
      "Finished sampling from posterior distribution ...\n",
      "\n",
      "Acceptance rate:  0.633\n"
     ]
    }
   ],
   "source": [
    "prior_k = HalfNormal(sigma=1)\n",
    "\n",
    "process_k = neuq.process.Process(\n",
    "    surrogate=neuq.surrogates.Identity(),\n",
    "    prior=prior_k,\n",
    ")\n",
    "# use burnin steps to tune to good acceptance rate, with random seed fixed\n",
    "samples, model = build_model(process_u, process_k, burnin_steps=2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "45bb1baa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEaCAYAAAAPGBBTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABU9ElEQVR4nO3dd3zT1frA8c9JmrRJOtl7yl4FykZkOFjXPcCJqLj3REUUvY6fePW6t6JeEecVryAiiiIOqOy9R9mjK22TZpzfHwk1QHez2j7v16uvNvmefL+naZonZz1Haa0RQgghSmKIdAWEEEJENwkUQgghSiWBQgghRKkkUAghhCiVBAohhBClkkAhhBCiVBIohBBClEoChRA1jFKqr1Lqd6XUL0qpmUopU6TrJKo3CRRC1Dy7geFa6yHADuCcyFZHVHcSKISoAKXUbqVUT//P9ZVS85VSmUqpd8v5+KeUUneEso5a631a6wL/zULA67/2EqVUl1BeW9RMEihEraGUSlFKaaVU6xPuf1Ep9VZ5Hg80Atb575oMbNZapwB3+8/9+wmPeV0p9bz/5/rAlcAbAcfrKKW+UkrlKaV2KqUuLefvcotSKl0p5VRKvV9CmZbAmcA3/rumA9PKc34hAkmgELVJKmDH1x0TqBuwohyP7wZs0lo7/bdPBz4LOPd+oLNSqlHAY3oGnHsCMCfg0z7AK/g+9TcELgNeK+en/r3AE0CxLRmlVCLwITBBa+3y3z0bGHZC/YQokwQKUZukAuv0yZkwuwDLAZRSDZRSs5VSB5RSuUqpb/xvugDdgdVKKbNSKhtf4PhGKbXaf+50YD7+MQGllNFfZrn/8aOAn49dVCllAy4Apmit7VrrX/G9mV9R1i+itf5Sa/1f4MiJx5RSMcAnwGNa640Bj3EAfwFnlXV+IQJJoBC1SU9gTeAdSqmGQF1glf+uROAloAXQEqgHXO8/1g1YpbUuBAYAB7XW8VrrbvzdcvgvcK6/fEd8/2PrAx5f9MYNtAfcWutNAfetxBe4qmI80A+YopRaqJS6JODYeqBHFc8vapmYSFdAiDBKBU5RSgXOAjICW7TWdgCt9RZgi/+YUyk1H0jx3+4O/C/gXCtPOPfXwI/A60qpBP99awO6fpKB3IDHxAM5J9QxG0io8G8WQGv9Ib5up+LkAo2rcn5R+0igELWCUioW6ARcDCwNOHQT0Dag3EXAHUA7wAxYgUlKKQV05e+WRyr+QBFw7hVa60yl1BJ83UyB4xMAmRwfBOz4WjCBEjk+mARbApAVwvOLGki6nkRt0RXf632+1jrj2Be+gHBsfGI48Ay+QNEEX7fTQXxv9q0Bj9Z6p/98Pfi7RdEVyAe2+W//F1/3U0/+Hp8AX5BpH3B7ExCjlGoXcF8PYG2VftPSdeL4lpAQZZJAIWqLnvgGsu0n3N+Hvz/198C3WG0lvu6md4EG+KbDduf48Y3AQNET39jFsUHy2cBoTm5RzAFOO3ZDa50HfAlMU0rZlFKD8A2EF3UbKaXeL276q1IqRikVh6/rzKiUivMPYpfIX743vgF3IcpNAoWoLVI5vsvp2LqGlvz9Zv4fwAQcxTcWsRlfcCnEP5Dtf1wjfIFkQ8C5j50DrfUOfFNwkzn+0/sHwGillCXgvpsAC76Wy0zgRq11YIuiObC4mN/nYaAAeAC43P/zwyX98n7/ABZqrfeWUU6I4yjZM1uI8FFKPYlvttQL5ShrxhdougcMiFfl2n8C12it15RZWIgAEiiEEEKUSrqehBBClCpigUIp1Vwp9ZNSap1Saq1S6vZiyih/Hp4tSqlVSqlekairEELUZpFcR+EG7tZaL/MvTvpLKTVfa70uoMwofNMX2+Fbafqa/7sQQogwiViLwp8KeZn/51x8qQWanlDsHOAD7fMHkKyUklWlQggRRlGxMlsp1QrfnPM/TzjUFN+89mMy/PftK+Yck4BJADabrXfHjh1DUlchhKiJ/vrrr8Na6/rFHYt4oFBKxQNfAHdorU/Me1NuWus3gTcB0tLSdHp6epBqKIQQNZ9SamdJxyI668m/l+8XwH+01l8WU2QPvgVHxzTz3yeEECJMIjnrSQHvAOu11v8qodhs4Er/7Kf+QLbW+qRuJyGEEKETya6nQfg2aFmtlFrhv+9BfPsAoLV+HV9unNH40j7nA1eHv5pCCFG7RSxQ+HfzUmWU0cDN4amRECIUXC4XGRkZOByOSFdFAHFxcTRr1gyTyVTux0R8MFsIUbNlZGSQkJBAq1at8PU4i0jRWnPkyBEyMjJo3bp1uR8nKTyEECHlcDioW7euBIkooJSibt26FW7dSaAQQoScBInoUZm/hQQKIYQQpZJAIYSo0Y4cOUJqaiqpqak0atSIpk2bFt0uLCws9bHp6encdtttZV5j4MCBwarucYYOHUpZi4dfeOEF8vPzQ3L9Y2QwWwgRVTweD3PnzmX58uX07NmTUaNGYTQaK32+unXrsmLFCgAeffRR4uPjueeee4qOu91uYmKKfytMS0sjLS2tzGv89ttvla5fVb3wwgtcfvnlWK3WkF1DWhRCiKjh8Xg466yzGD9+PFOnTmX8+PGcddZZeDyeoF5nwoQJ3HDDDfTr14/77ruPJUuWMGDAAHr27MnAgQPZuHEjAAsXLmTs2LGAL8hMnDiRoUOH0qZNG1588cWi88XHxxeVHzp0KBdeeCEdO3bksssu49jmcHPmzKFjx4707t2b2267rei8gQoKChg3bhydOnXivPPOo6CgoOjYjTfeSFpaGl26dGHq1KkAvPjii+zdu5dhw4YxbNiwEstVlbQohBBRY+7cufz555/Y7XYA7HY7f/75J3Pnzi32jbUqMjIy+O233zAajeTk5LBo0SJiYmL44YcfePDBB/niiy9OesyGDRv46aefyM3NpUOHDtx4440nrUdYvnw5a9eupUmTJgwaNIjFixeTlpbG9ddfzy+//ELr1q0ZP358sXV67bXXsFqtrF+/nlWrVtGr199b8Pzzn/+kTp06eDweRowYwapVq7jtttv417/+xU8//US9evVKLNe9e/cqPVfSohBCRI3ly5eTl5d33H15eXlFXUfBdNFFFxV1aWVnZ3PRRRfRtWtX7rzzTtauXVvsY8aMGUNsbCz16tWjQYMGHDhw4KQyffv2pVmzZhgMBlJTU9mxYwcbNmygTZs2RWsXSgoUv/zyC5dffjkA3bt3P+4N/tNPP6VXr1707NmTtWvXsm7dumLPUd5yFSGBQggRNXr27InNZjvuPpvNRmpqatCvFXidKVOmMGzYMNasWcM333xT4jqD2NjYop+NRiNut7tSZSpq+/btTJ8+nQULFrBq1SrGjBlTbB3LW66iJFAIIaLGqFGj6NevH/Hx8SiliI+Pp1+/fowaNSqk183OzqZpU9++ae+//37Qz9+hQwe2bdvGjh07AJg1a1ax5YYMGcLHH38MwJo1a1i1ahUAOTk52Gw2kpKSOHDgAHPnzi16TEJCArm5uWWWqwoZoxBCRA2j0ci8efOYO3cuK1asIDU1tcqznsrjvvvu46qrruKJJ55gzJgxQT+/xWLh1VdfZeTIkdhsNvr06VNsuRtvvJGrr76aTp060alTJ3r37g1Ajx496NmzJx07dqR58+YMGjSo6DGTJk1i5MiRNGnShJ9++qnEclWhjo3I1ySycZEQ0WP9+vV06tQp0tWIOLvdTnx8PFprbr75Ztq1a8edd94ZkboU9zdRSv2ltS52LrB0PQkhRBi89dZbpKam0qVLF7Kzs7n++usjXaVyk64nIYQIgzvvvDNiLYiqkhaFEEKIUkmgEEIIUSoJFEIIIUoV0UChlHpXKXVQKbWmhONDlVLZSqkV/q9Hwl1HIYSo7SLdongfGFlGmUVa61T/17RQVqbA48FbA6cLC1HbGY3GotTiqampPP3000E794oVK5gzZ07QzheNIjrrSWv9i1KqVSTrEMjp9ZLj8dDQbI50VYQQQWSxWEKSLwp8gSI9PZ3Ro0eH5PzRINItivIYoJRaqZSaq5TqEuqL7Q5CXhQhRPTLzs6mQ4cORSnFx48fz1tvvQWUnKp76dKlDBw4kB49etC3b1+ys7N55JFHmDVrFqmpqSWm5qjuon0dxTKgpdbarpQaDfwXaFdcQaXUJGASQIsWLSp9wQMuF4VeL2ZDdYihQlQvd9xxR9A/2aempvLCCy+UWqagoOC4xIKTJ0/mkksu4eWXX2bChAncfvvtZGZmct111wHFp+ru2LEjl1xyCbNmzaJPnz7k5ORgtVqZNm0a6enpvPzyy0H9vaJJVAcKrXVOwM9zlFKvKqXqaa0PF1P2TeBN8KXwqOw1vVqzx+mktcVS2VMIIaJMSV1PZ5xxBp999hk333wzK1euLLr/008/5c0338TtdrNv3z7WrVuHUorGjRsX5WlKTEwMV/UjLqoDhVKqEXBAa62VUn3xdZUdCfV1d0mgECIkyvrkH25er5f169djtVrJzMykWbNmRam6ly5dSkpKChMmTAhKqu7qLNLTY2cCvwMdlFIZSqlrlFI3KKVu8Be5EFijlFoJvAiM02HIYpjjdpMThBzyQojo9vzzz9OpUyc+/vhjrr76alwuV4mpujt06MC+fftYunQpALm5ubjd7uPSfNdUkZ71VPw2T38ffxmISMffbqeTLiVsuC6EqF5OHKMYOXIkV199NW+//TZLliwhISGBIUOG8MQTT/DYY48Vm6rbbDYza9Ysbr31VgoKCrBYLPzwww8MGzaMp59+mtTU1KKxj5pG0owHyHK5WJSdDYDZYOCMlBQMSgW7ekLUKpJmPPpImvEgKfR62VdYGOlqCCFExEmgKMXOWj6AJYQQIIGiVEdcLuwyqC2EqOUkUJRhp9MZ6SoIIURESaAow26HA08NHPAXQojykkBRBpfW7JNWhRCiFpNAUQ7S/SRE9aaU4vLLLy+67Xa7qV+/PmPHjo1graoPCRTlcNTlkpXaQlRjNpuNNWvWUFBQAMD8+fNp2rRphGtVfUigKKftMlVWiGpt9OjRfPvttwDMnDmT8eP/TgyRl5fHxIkT6du3Lz179uTrr78GYMeOHZx66qn06tWLXr168dtvvwGwcOFChg4dyoUXXkjHjh257LLLqImLl4+RHBXltMfppLPViknSjwtRNUOHnnzfxRfDTTdBfj4UtwHQhAm+r8OH4cILjz+2cGG5Ljtu3DimTZvG2LFjWbVqFRMnTmTRokWAL6348OHDeffdd8nKyqJv376cfvrpNGjQgPnz5xMXF8fmzZsZP348x7I+LF++nLVr19KkSRMGDRrE4sWLGTx4cLmfhupEAkU5ebRml9NJW8kqK0S11L17d3bs2MHMmTNP2o3u+++/Z/bs2UyfPh0Ah8PBrl27aNKkCbfccgsrVqzAaDSyadOmosf07duXZs2aAb49MXbs2CGBQsAOh4M2cXEoyf8kROWV1gKwWks/Xq9euVsQxTn77LO55557WLhwIUeO/L1jgdaaL774gg4dOhxX/tFHH6Vhw4asXLkSr9dLXFxc0bHY2Niin41GI+4aPI4p/SgVkO/xcEDyPwlRbU2cOJGpU6fSrVu34+4/66yzeOmll4rGGZYvXw74tktt3LgxBoOBDz/8EI/HE/Y6RwMJFBW0Qwa1hai2mjVrxm233XbS/VOmTMHlctG9e3e6dOnClClTALjpppuYMWMGPXr0YMOGDdhstnBXOSpImvEAgWnGSzM0OZkE2atCiHKRNOPRR9KMh8E2aVUIIWoRCRSVkOF04vR6I10NIYQIi0jvmf2uUuqgUmpNCceVUupFpdQWpdQqpVSvcNexOF6tZaxCCFFrRLpF8T4wspTjo4B2/q9JwGthqFO57JCsskKIWiKigUJr/QtwtJQi5wAfaJ8/gGSlVOPw1K50hV4vGZIsUAhRC0S6RVGWpsDugNsZ/vuiwraCghqd30UIIaAGrcxWSk3C1z1FixYtwnJNu8fD/sJCGges0BRClO6bw4eDer5/1KtX6vHdu3dz5ZVXcuDAAZRSTJo0idtvv71c516xYgV79+49KeXHMa1atSI9PZ16ZdShuov2FsUeoHnA7Wb++06itX5Ta52mtU6rX79+WCoHsElaFUJEtZiYGJ577jnWrVvHH3/8wSuvvMK6devK9dgVK1YwZ86cENcw+kV7oJgNXOmf/dQfyNZa74t0pQLluN3sl7QeQkStxo0b06uXb8JkQkICnTp1Ys+ekz9vfvbZZ3Tt2pUePXowZMgQCgsLeeSRR5g1axapqanMmjWLI0eOcOaZZ9KlSxeuvfbaog+JeXl5jBkzhh49etC1a1dmzZoV1t8x1CLa9aSUmgkMBeoppTKAqYAJQGv9OjAHGA1sAfKBqyNT09JtKiigkdksyQKFiHI7duxg+fLl9OvX76Rj06ZNY968eTRt2pSsrCzMZjPTpk0jPT2dl19+GYDbbruNwYMH88gjj/Dtt9/yzjvvAPDdd9/RpEmTov0ussuR4aE6iWig0FqPL+O4Bm4OU3Uq7VirQsYqhIhedrudCy64gBdeeIHExMSTjg8aNIgJEyZw8cUXc/755xd7jl9++YUvv/wSgDFjxpCSkgJAt27duPvuu7n//vsZO3Ysp556auh+kQiI9q6naiOSYxVaa1xeL06vF4fHg1fGTIQ4jsvl4oILLuCyyy4rMQi8/vrrPPHEE+zevZvevXsfl4a8LO3bt2fZsmV069aNhx9+mGnTpgWr6lGhxsx6irQct5sMp5PmAfnqQ63A42GX08kuhwNHQEoRs8FAE7OZprGx1DGZwlYfIaKR1pprrrmGTp06cdddd5VYbuvWrfTr149+/foxd+5cdu/eTUJCArm5uUVlhgwZwscff8zDDz/M3LlzyczMBGDv3r3UqVOHyy+/nOTkZN5+++2Q/17hJIEiiNbm5dHAbCY2xNulFnq9rM3LK3HBX6HXyw6Hgx0OB3VNJrrZbJLtVkSNsqazBtvixYv58MMP6datG6mpqQA8+eSTJ015vffee9m8eTNaa0aMGEGPHj1o0aIFTz/9NKmpqUyePJmpU6cyfvx4unTpwsCBA4um4q9evZp7770Xg8GAyWTitdeiJolEUEia8QDlTTNemqaxsfRKSKjSOUqz1+lkdV4ehRVISqiAthYL7a1WjDLgLsJM0oyHj9a6XJNqJM14hO1xOkOyC16h10t6Tg5/5eZWKEgAaGBLQQG/ZmeTV0t36BKipvNoTahyWkugCIHVdjvuIKYhP1RYyM9ZWeyrYgDKcbv5JSuLvZKjSogaJ5RbH0igCIECr5f03NwqZ5d1eb2ssdv5IyfnuMHqqnBrzV+5uazPywvK+YQoj5rYxR1NvFpTWM7nuDJ/CwkUIXLI5WJJTk6lg0WGw8FPWVlsD9G+F1sKCliRmytTaUXIxcXFceTIEQkWIeT0eqEcz6/WmiNHjhBXwdmZMhUmhA67XPyZk0PfhARiyjETyqs1e51OtjscZLndIa/fbqcTp9akJSTIILcImWbNmpGRkcGhQ4ciXZUaSWuNw+tFA3EGA4Yy/pfj4uJo1qxZha4hgSLEjrhc/JiVRYvYWFrGxWExGo877vJ6Oep2c6iwkAynE1eYP3UdLCxkaU4OfRMTy3yBCVEZJpOJ1q1bR7oaNdbWggK2+ruShyUnEx+CqfASKMLA6fWyuaCAzQUFpMTEoJRCa41ba3KjYBbSIZeLlXY7PUM4rVcIEXxerdlWUBDy60igCLPMMHQpVUaG04nFYKCjzRbpqohqpNDrZX9hITluNwVeLwVeL7EGA4lGI0kxMTQwmcrV7SoqZ7fTGbSJLqWRQCGKbC4owGo00iKMaUhE9bTP6WSnw8Fhl4viOksP+r+bDQZOsVhoFRcn42BB5tWazfn5YbmWBApxnDV5edSJiQlJP6eo/nLcbtbk5XHE5SpX+UKvl3V5eWwtKKC7zUYjybAcNBlOJwVhaE2ATI8VJ/BozV92u0ybFSfZkJfHL1lZ5Q4SgZxeL0tzc1mblyevrSDwas3mMIxNHCOBQpwkx+1mfZiatCL6aa1ZkZvL5oKCYruZKmJbQQGLs7MpiIJJHNVZhtNJfhifQwkUoljbCgo4JFu81nperUnPzWV3ENO+ZLndLMrOJrMSLRMR/tYESKAQpVhpt1c5DYmo3pbb7SHZE97p9fJbTg4ZIco8UJPtdDjC2poACRSiFAVeLxulC6rW2uVwhDSBpFdrltvtMm5RAS6vl01hbk1AhAOFUmqkUmqjUmqLUuqBYo5PUEodUkqt8H9dG4l61mbbCgrIidK1HyJ07P7ZTeGwraCA33NycMi4RZk2FxRUeJuBYIhYoFBKGYFXgFFAZ2C8UqpzMUVnaa1T/V81a3/BakDj64KShG61h9c/8y2c3Y5HXS5+yc7moIyLlSjf4wlZktCyRLJF0RfYorXeprUuBD4BzolgfUQJstxudkhfcq2xMT8/Iq1Ip9fLnzk5LK/E5lyVke/xkO12k+12k+N2R32LZkN+fsS66CK5qqopsDvgdgbQr5hyFyilhgCbgDu11ruLKYNSahIwCSjax1YEz4b8fBqbzcSdkNRQ1CwFHg/bIvyhIMPp5JDLRWerlaaxseXa2rM8Cr3eonNnud3FBiOTUsQbjSTHxFDPZKKuyYQpClKQHHG52BPBDceiffntN8BMrbVTKXU9MAMYXlxBrfWbwJvg2zM7fFWsHdxasyYvj7TExEhXRYRQJD+1BnJ6vSy329nqcNDJaqWB2Vzpcx1xudjpcLCvsLDM382lNZluN5luN9sdDhRQz2SiZVwcDc3miGRY9mjNCrs97NcNFMlAsQdoHnC7mf++IlrrIwE33wb+Lwz1EiXYV1jIfqez2qRhcHg87PDv7ZHn9eLwekmOiaGR2UxjsxmrtI6Ok+N2kxFl2+TmuN38mZNDSkwMLePiaBIbW66cUV6t2VdYyNaCArKr0I2m8WVXPuRyEWsw0NZioXVcXFgDxvq8vLBPhz1RJAPFUqCdUqo1vgAxDrg0sIBSqrHWep//5tnA+vBWUZxodV4e9aI8I2iBx8OWggJ2OZ0nfYI86nJx1OViXV4e7a1W2lssQevaqO6ieTV+pttNpn8qbZPYWOrExJBiMmHzB/tjafsPu1zsLyzkoMsV9HEOpz9v1S6Hg242G/Wq0MopryMuV8QGsANFLFBord1KqVuAeYAReFdrvVYpNQ1I11rPBm5TSp0NuIGjwIRI1Vf4OLxeNuTn0zU+PtJVOYnWmu0OBxvy88s1Y2dTfj5HXC56xcfX+rGXw4WF1WLGkUtrdjoc7PTfNiqFhrB2l9k9Hn7PyaFVXBxdbLaQtS7cXm/Eu5yOiegYhdZ6DjDnhPseCfh5MjA53PUSpdvucFDfbKZhGD5RlVeO280Ku73C3QxH/NMyByQmklCLM+ZujMAirmCIZOaAHQ4Hdo+HtISEoA94a61ZZrdHvMvpmNr7n1FD5eXkkLF5M0f37yfzwAGUwUDzdu1o3qEDSfXqBe06K+x2TktKivgncY/WbMrPZ2sVEtY5vV7+yMlhUFJSrRy3yPR3x4mKO+xysSg7m36JiUXdYMGwLj+fA1HUwpNAUQM48/NZMm8eP3/5JcsWLMBdwj99o1atOOvKKznjsstIrFOnStcs9M9K6Z+YGLE+/iP+LVzzSvnU5Xa5OLp/P4f37iXr4EFadu5M07ZtTyrn8Hr5PSeHQYmJEQ9+4balmrYmokWex8Nv2dkMTEoKSrDYXlAQlu1NK0LVxBW3aWlpOj09vcKPy/J/OqguPG438z74gJnPPkv24cPUadSIU889l26DB1O3cWPqNGqEu7CQ3Zs2sXvTJpZ89x2rFy/GFBvLsIsv5sqHH65ywOhgtdLeag3Sb1Q+bq+X9fn5pS4C3L9zJ/999VV+mDmTwhP+6Rq3bk2fM89k9MSJNGnT5rhj8UYjg5OSomLufDjkeTz8mJkZ6WrUCHEGQ5WDxV6nk2W5uZVuHQ9LTq70pmNKqb+01mnFHpNA8bfqFCiW/fgj70yZwu5Nm+g6cCCX3HMPXQcOxFjGi3Tn+vV8++67zP/oI+JTUrjh6acZdPbZVapLj/j4sG2fesTl4q/sbH79/nu2rl5N227d6DViRNHvnXnwIO8/8gj6q6+4GDjXYuGdW2/F07Mng37/ncZz5vC50cgrW7eSGxPD9c88w/BLLjmuVdTQbKZPQkKtmA21ym5nZxTMqqkpqhIstuTnV3nmmQSKCqjJgaLQ4eD9xx7jf2+/TZM2bZjw6KP0Gzmywm9q29eu5aXbb2fLypUMHDuW2158EWtCQqXr1SshgaYhXF+htWZTQQHrc3OZetFFbFq2DGd+PrFWK+179eKxzz5j3e+/k3HllTycm0sTwB0Xx6Hhw1nz9NM4Gjem3sKF9Lj7bqy7duGJieHj+vWZtG8ffc8/n5umTz/u949ESyncnF4vP2RmRsUCu5ok1mCgX2IiSeV8w9ZaszovLygBWwJFBdTUQLF70yaenTSJHWvXcvb113PVlCmYqvDm7HG7+eqVV/joqado3r49D3/0EQ0rmf5EAb0TEmgcgmDh8Hj4y27nqMvF0u+/59lJk3AEZDaNtVoZMGYMuz7/nPVac6RLF/bccQcHzzgDj812/Mm0JmnFClp+8AEtP/qI1a1akbp7N+169uSJL78k1mIpKtovMbFKK4Kj3Ya8vLBvgFNbxChFWkIC9ct4/WS73ay228kMUm4tCRQVUBMDxR9z5/KvG27AbLFw+4sv0ufMM0surDW2rVtpOH8+1u3bUcD6KVNwJyQQk5uLx2JBB7yYVvz8M89MnEiM2cxDH3xAxz59KlVHBXS0WmkbxEVsBwoLWWG3Fy2e+uS555j5zDPHZbNNBVYAp11wAQ9ccQW/5+ayZe3ak7qlTlT/xx/RRiPf2O08ffXV9Bs9mvvfeaeovEkpTktOxlIDB7c9WvNDZmZEUlbXFgrobLPRPDb2pDEvl389UrCTbUqgqICaFCi01nz58st88PjjnJKayoMzZlC3ceMSyzf+73/p+OSTxG/fDkBhSgraYGDhokUU1q9P50ceofnHH3Ng5Eg23nsvBS1bApCxeTPTLr2UI/v2Mfm990g744xK17mR2UxqfHyVBoS9WrM+P/+k2R8ntiimANOAyQMG0PfLL3n04ouL7ZYqdexGazzjx/P5ggXkXn891z7xRNGhBmYz/WpgfqtdDgcro2QxV01nUIoGJhMNzWbsHg9HXC6y3e4q7z9enFAFitoxtaOachUW8uLttzNj2jQGnX02T379dfFBwuslJjcXgNwuXchv3ZpVzzzDD8uWMW/TJr7fsIHC+vUBODh8OPtHjaLJ7NkMGziQTo8+Skx2Ns3atWP6d9/RokMHnrzqKv6YO7fS9d5fWMgv2dnsr2TeoBy3m8XZ2cVOEew1YgTte/XCZrHwDr4gMSs+nv6ffsryH39k07JlOPLy0FrjyMtj07JlLFuwoNTrKbebwU4nM5RCvfEGc959t+jYwcJCdtfAwd5om35Zk3m1Zn9hISvtdrYWFJAVoiARShIoolS+3c7jl13GgpkzueTuu7nnzTeP6z8/xnzwIAPOO48et90GgL1dO/6cNYudEydS0Lz5SeUPDx3Kypde4sc//2TP+efT9tVX6faAb3PBxLp1efzLL2nTrRvPTJzI4tmzK19/j4elubkszs4mq5yLudxeL2vz8vglK4usEvpsjUYjj3z0EbPi4pgIfJSaSuzmzRji4ti6ejXOE2aNOPPz2bZ6danX1SYTSz78kOyePflUKbY9/DB7t24tOr42Ly/q9yqoiEOFheTWoN9H+DScNw9C9HeVQBGFju7fz4Nnn82qRYu47d//5rIHHsBQTDdO4urVDDnzTFKWL+fgiBFQgW5ER+PGrHzpJX754QfWPPkkADG5ucQnJjLt889p17Mnz06axG/ffFO138XfnfdrVhbbCwpwntAnrrUmy+ViQ14eP2Vlsa2MFdZer5ffrrySMZmZfDFyJEnz52P0N7XbdutG7AkzlWKtVtp061ZmPT3x8Sz55BPy2rblP243H998M15/XV1asypM24KGQzQkmRPBY3A46H7nnfS9/HJiPvggJNeQMYoA0TBGsXvTJh4bP56cw4e5/5136H366cWWa/Ttt/S86SZcycks/eADsnv0qNJ1VWEhA885B0eTJix/+WXyPB6mXnwxm5cv58EZM0ofPK/IdfBNH4xRCpNSFPjTf5eH1pq3HnyQ/739Ng9OmED/Z5897rjH4zlu6qzFaqVX37588M03uJQi3+PhoMt1UrAKFL9pE/1HjOByh4OCq68mpWHDokHx/snJ1SbFeklkgV3NYzp6lCFnnEHGBRfQ5KmniK/ka1QGs8sp0oFizW+/8eRVVxFjMjHl449pl5pabDmj3c7wfv0oaN6cpTNm4GzYsOoX15q2r7xCp2nTODx4MEs/+IAcr5cpF1zAjnXrmPKf/9Bz6NCqX6cKDo0fz4c//EDSDTcwcdq0k2ZWGZWirdnMxoULWbVyJampqYwaNeq4gWytNYdcLjKczhJ3DFNHjjC+f3/sWVkAxNlstO/Vi+lffsmIunUjsnlNsKyx24PWotAOB1mrV7Px0CFirVZaxseT2KkTMSdOSRbBpzWNv/mGA2edhTc2FqPdjic+XmY9VUR1DBS/fPUVL9xyC41atuSRmTNp5J+NVBLb1q04GjU6eZ1AFTX79FNSb72VzN69+fOTT8j0eHjovPPYu20bj86cSddBg4J6vfLaeeed3PTRR/zUsiV5S5eeFCTqm0x0j4+vUFK/TJeLFXY79hP6dZd+/z3PXHstowoK2AqswRcs7n3zTSaedx6tihkrqg5cXi/zMzOrlHHVaLcT889/0ubjj2mdn8+fwLFXxGqgE7A6IYGtY8dimjwZVcoMPVE5cXv30v2ee2g4fz6rnn2WnRMmFB2TWU81lNfrZeazzzJ90iTa9+rFM99+W2KQaDB/Ph2eegq0Jq9t26AHCYCMiy/mr7feInn5clJvu42ElBSmffYZDZs3Z9qll7Luzz+Dfs2ybPzXv5jw0UdsSEoi/8cfTwoS9Uwm+iYmVjjza4rJxGnJybQ74Y1/6+rVmAoKeA3f3rsG/h4U31RQgLuarj3Y4XBUKUikPPQQp7Zrx5lvv02my8WnqamsmzCB6fPm8dTs2Sy/8kq+7tOHOIeDC2bOZGT37hScdx5OmWEVFDHZ2XR84gmG9+9PvV9/Zc0TT7DzyivDc+2wXEUUK99u54Wbb+aPOXMYfskl3DR9OuYSciYlp6eTNnEiue3bs+X22/GEML3EvrPPZonNRp4/y2py/fo8/sUXPHTuuTw2bhyPffYZHdOK/eARdOtmzuTSp57CbjazfcECOGFNQ1JMDH0SEirdHWRQio42G2aDgbX+Aeu23brhsdm4PS+PT4FrgRkWC226dcPp9bLV4aBDNUvv4fVv6lRZ382Ywfa33uKamBjWX3stHR95hHiLhXigqONzwAAANno8LProIxJffpnXf/2VdYMHc8uUKfQ+/XQ8UbjhVXlprcnYvJmN6els/Osvdm3YgD0ri9zMTAodDmxJScQnJ5Ncrx4tOnWiTbdunNKjB83atcPr9bJswYJi85OVV89bbqHRd9+Rcf75bJw8mfxWrULzixZDup4ChLPradeGDTxz7bXs2bKFiY89xj8mTSpxNbNl1y5OHTkSt83Gr3PnUhjEfSXK5PXS+u232XX55RzMzmby2WeTc+QIj376aciDxZ/ffUfsVVdxF/DrN99Q2Lfvccet/kyvsUHK9LoxP59N+fl/D4r/9Rdz8/PpCoxo357HfvkFo9GIUSmGJydXq3Tkux2OSu2W1uSLL1g7cyZX/fwzaaefzl2vv058UlK5H7968WJev/9+7ti4kXEWC6vffZfCEiZoRCNDQQG5v/7KV+npLPz8c/65axcDgBSlsBqNFJhMbG7QgOlnnUVedjaD1q4lKyeHn/ftY63LhQdIadgQtMaenY3L6Swa8yptIajRbqfe4sU0++QT1j36KAUtW5Kwdi3K6yWnlFl8MkZRAdEcKLxeL7PfeIMP//lPrAkJ3PPmm/Q49dQSy8fk5DB49Ghi9+9n8dy52Nu1C2n9TpSydCmDxozh0NChLP3wQw4cPsxD555L1sGDTJ4xI2QD3H9+9x3PTJxI265dee6pp/D27n3ccYNSDElKCvqudMcGez0eD8sWLMA5fz5Pvv8+HyYlkbxhw99TcS0WOlejQduFmZkVXjvR8o036P7wwywEnrn0Uq5/7rmi378i3C4X6yZPZvyMGbQEfr/kErJeeOG4NDLRxFBQQIP587G8+y6tf/+dLK+XFgYDPYYM4RmDgdZaY2rcGG9sLDF5eRQ0bcrGBx8EYETv3lh37QLAbTazr0EDPleKu3bvBqA3vj2ds0wmbnruOQafdRZesxlPfDymrCzavPYadX//nZT0dAwuF866dVnx8sscLGdwrZGBQik1Evg3vj2z39ZaP33C8VjgA3zP7xHgEq31jrLOW+lAkZ3NrrfeYucVV0AIPi3u2bqVV++5h9W//krfkSO55V//Itm/Yrok9X/8kbSrr2bJRx9xpJSAEkrN//MfUu+4g/2jRpH+zjscPXqUqZdcQsamTdz9+utVTlN+ot+++QbvtdeS3qkT1379dbGfYDvbbLQNwaCy1prfcnKO2/HNeu21zPj6a5Kef54zL78c8CV9Oz0lpVrsW3GwsJA/c3Iq9Jhms2bR85Zb+Ar44q67uPiBB6qcv2vXb79R77LLuMRuZ1OTJmyfM4fCpk2rdM5ga/Lll3S6+26sdjuHgDkWC1nnnEPSgw+WmjrnGGN+PrZt20hct46EdetIXLuWX10uLli8GAUUcnJ//4/DhrHvtddIiY1lVNu2ZHftyuHTTuPQaadxtH9/vBWY7lrjAoVSyghsAs4AMoClwHit9bqAMjcB3bXWNyilxgHnaa0vKevclQ0U+a+/jvXGGznSvz/LX365KA9SVWUePMgn06cz74MPiI2L45onnuCMyy4r9z+e+fDh8HY3FaPVO+/Q7YEH2HPOOSx//XVy8/J4/NJL2ZiezlVTp3LujTdW+Y1Ea81nL7xAnSef5N/AmptuYvtjj51Uro7JxMAQ7qxnd7v5OTu7KP221pr7Ro/myL59vLFkCSZ/RtDqkor89+xsDldgq9OGc+fSe8IEFnq9vDJhAlefsF6lKvLtdjaOG8e1f/7Jrf37c8EHH5CQkhK081eGJSMDr8nE5qwslt12G2OXLePj+vVpdO+9DLv00iplaIa/85MV5uVxFlAfaBwTQ9N27di/fz8/Z2byu8FAxz596H/66XQbPpxWXbpUeAwDIhgolFLpwEp8s99WA6u01ocqVZPjzzsAeFRrfZb/9mQArfVTAWXm+cv8rpSKAfYD9XUZla5soLho3DgGHTrE9b/9hlEpVk2fzr6LL67wecD35rJt1Sp+mDmTHz/5hEKnk7OuvJJL7r6blAYNynx8s1mz8MTFse+ccyp1/VBo8+qrdHzySRZ/+y3ZPXrgzM/n+Ztv5rf//Y/+Y8Zw27//XaH+60CFDgcv3XknHT7/nDeAPSNHsuz9909q2Rn9GV2DuT9xcU7cRGbFd99x8IoraHfrrbR+5BHAl1329JQUYqK4VVGZ7lTLY48R/8or3NGvH/d/+SUxJlNQ66S1ZsE77/DqI49Qp1Ej3rj+eoyTJkGY16eYDx+m3Qsv0PK99/i+bVvGbtqExWZj3L33MnrixKIPBFV14kLQwGSVBoOBratWseS77/jzu+/YvmYNAJb4eDqmpdG6Wzcat25NkzZtSK5fnzirlVirFe31kpedTV5uLkf372ff9u3s3baNBkoxIyBXWUVUNVA0Abr7v9KAMcBhrXWVPm4rpS4ERmqtr/XfvgLop7W+JaDMGn+ZDP/trf4yh4s53yRgEkCLFi1679y5s0L1cblcDBw8mFWrVtHI4eBDYAgwtXlzlowYQatOnWjSti1NTzmFuo0bn5RSQ2tN5oEDbFq2jE3LlpH+ww/sWLsWU2wsA//xD8bdc0+xezUXJ2XpUgacey5HBw7kj08/Dfs/UGksO3f+3dLSGg3MfuMN3n/sMeo3a8ZtL7xQ4bUWa377jdfvv58RGzbwLrD/jDNIf/99dDH/qN1strCsY/Bqza/Z2WT7c04ph4Oebdpw0GBg544dRX31XWw22kTxuor0nBz2FRaWu3zmwYPcOXw4sWYzz86fT2LduiGr28a//iL90kv55OhRVnTtyt6vvsKTnByy6x2jCgtp/dZbtJ8+HWNeHh/HxvKQw0Hnq67isgceICkErfdjY17bVq+mTSmzng7v3cvaP/5g/R9/sG7JEjI2b8Zdzr+fLSmJHt268esvv1SqtR3UrielVCfgQq314xWuyfHnCWqgCFSVweyFR46wZ8sWdixfzmlvvcUrZjOLNm8mL6CP12AwYE1MxJaUhNFoJC8nh7ycnKI/qDEmhnY9ezL0oosYct55xFfgxR+3Zw+nnnkmHquVRd9/jyvCzfKSNJs1i6aff076e+/hiY9n/ZIlTL/+eg5lZJB2+ulc8fDDtO7SpdRzHNm3jw+eeIKfPv2URs2a8afFQlyzZiz94AO8xUwTTomJYVBSUti2KM12u1mUlVWUe8p5991c+MEHzLj1Vur4WxVxBgMjUlKicrV2hdJ1aE23u+/mlfXr+ffq1UyfN6/Mv18wZB86xP6zz+b2LVs4aLOx9pNPcPbvH9Jrdn7kEdq+9hq/1KnD9UeP4u7RgxuffZZ2PXuG9LqV4fF4OLxnD3u3bSM3MxNHXh7O/HwMRqPvPSgxkaR69WjSpg0JKSkR7XpqqbXeecJ9H2qtr6hUbf4+R9R1PZXUTNceD43+9S8W9uzJjowMju7fjz0rC3t2Nl6PB1tSEraEBFIaNqR9r1607tq12EyvZTHa7QwaOxbrrl38Oncu9g4dKnyOcGk2cyY97ryTnG7dWPLRRzgbNsRZUMD/3n6bz//9b/JzcugyYADdhwyhx6mnUrdxYwodDpz5+axfsoTF33zDuj/+INZo5OLrruMfDzyAzePBGxODt5jnTgFDkpNJDPNMmZV2O7v86w+8TifdW7XCERPD5h07UP5PhOHcM7wiKrIfdosPPqDH3XdzP3BgyhQu8GcjDgev18uKe+5h4ocf0hD4+YYbcDxepc+hJ0lcswav2cyRZs34Ydo0dr//Pj/bbFw2eTKjrr66UuMB0SiSgeI3oAWwHd8YRRYwSmtdpfDrf+PfBIwA9uAbzL5Ua702oMzNQLeAwezztdZlDhoEO1Ak//UXg0ePJrN3b5Z8/DGuEDWPm3/0Ed3vuYclM2dyaNiwkFwjmBrOm0fv667DY7Gw6v/+r2g8xZ6Vxew332TpvHlsW72a4l5jLTp2ZNSwYUz9/XdikpL447PPSu1ia2Ox0CUC01ELPB5+zMoqGtg+cuedTPjoIz6/7TZip0wBIN5oZGhycthaOuVRkf2w4zdu5NQRI/jF6+WGTp34v3nzKjUNtqq2zJ1Lqxtu4N38fIw33MDlDzxAbBX/5pZdu2j3r3/R4uOPWd29O6cfOcKhjAyGXXwxE6ZOLdd4YXUS0VlPyvcf0BboBtQB5h3rDqoKpdRo4AV802Pf1Vr/Uyk1DUjXWs9WSsUBHwI98U0/Hqe13lbWeUOxjqLxN9/Q6/rryW3fnj8+/ZTCEL3A4jdswN6xY0jOHQrxmzeTevPNpCxfzs8LFpDTvftxx3OOHmXN4sXk5eRgtlgwx8bSukkTTl2wgLavvorB7WbFCy+w9/zzS7xGnMHAsOTkiA0ar83LK9rox+Vw4Grfnv+1acPFCxcWlembmEjDKNpfe31eHlvKkzrD62XQ2LEYV66ki9fL/QsW0Kpz59BXsAT5djszpk1j7nvvMT0piaEdO5L53HMUVLB1bdm5k3YvvEDzTz7BC7yflMS9R45Qr0cPJk6bRteBA0PzC0RYjZseG0qhWnBX/6efSJswAWeDBvzp37sgGNq8+ipHBg2qcqrwSFFuN/V//JGD/lTkbV98EXuHDhzt08fX+jIYfHtlKIVt61YGjRlD7JEj7Bszhg0PPVTmIsK0hAQaRzC9d6HXy4LMTNz+/5X/vvYa7z7yCM/Nn1+U4beOycSgSs74Cja3vzXhKsf/doPvv6ffZZcxASi87z7G33tvyOtXHqt+/RXrpElcf+gQscDqPn1wTJrE0VLSgMTk5uI1mfDGxdFy+nQ6TZ/Of+LjeTg7G0+rVlzx4IMMOuecYvd2qSkkUFRAKFdmJ//1F2lXX83yV1/lyODBla1ikXbPPUfHp59m+8SJrHnmmSqfL9IMDgfD+/XDsncvANpgwJWczNYbbmDLnXcSk5ND7+uuY8MDD5BdjsHDRmYzfaJgz+pj6T0A8nNzub9rV+5v25ZmP/zgC4TA4KQkUoI8lbQythUUFOWtKou7sJBZffrwa1wcLy5aFLQpocHg9XpZ+eGHNHjiCcZlZZEMzG7Tht+nTOGUpk057aOPUPHxWLdvx7phAwk7d/LOhRfyVlYWexctQhcUYO3Rg39MmsTgc8+Nqt8tVCRQVECoU3gYHI6iWTnxmzdXOq3GsSCx++KLWfHiiyFZDR4JyuWi/k8/Ydu+HfPRo5iPHuVo//7sueCCCp0nRimGJidjiYLn5cRP6ZvHjeOuBQv44fnnKfCv1m5sNpMW4aDm1ZofMzMpKEeGW/Phw3zx5Ze89dBDTPnPf4K2OVWwedxuVs6fT87777No6VKW5ubSEkgH4vENnq4H1gEzgaOtWpE6dCjDL76YDmlpUTV2FGoSKCogXLme6i9YQP9x49h16aWs+ec/y50Z0+B00u2ee2jxySc1LkgEU7StUdicn88Gf6vi0PbtDOnbFxo0YP2aNUUD8cNTUkK+GLA0e5xOluXmllmu/k8/kXbllYyKieFwWhqPffpptXhDdRUWsmXFCg7s2sXBXbs4un8/Vv8U0aR69eiYlkajMGZVLU1iTAynWCyYlcKgFFpr9hQWssfprFK699KEKlBEZ1auauLwqaey+Y47OOXFF6m7eDErn3/e1x1Vxj+c12zGfPQom+66i4333SdBohjJMTG0jrIpp63i4thaUIBLa+q3bs1nPXoweeVK9syZQ86YMYAvmKQmJESsjlvLMYCtXC66PPwwB8xmfs3N5dlidguMViazmU59+9LphEzC0cRmNNLRaqWx2Xzy3ilmM52tVnY5nWwpKKCwmuxtUnNHdcJAm81seOghFs+eDUox8PzzSQvYbeoYY14eTb78kj5XXEHC2rWgFEs//JCNkydLkCiGwrc2IdrevEwGw3EtHPPjj7MLaDp1qm+wHshwOsmrYJbWYDlcWFi0krw0LWfMIGHTJm6y2xk2YQItO3UKQ+1qh1ZxcQxNTqZJbGyJr1+TwUBbi4XhycmcYrFE5WLNE0mLIggy+/Xj54ULafbZZ3j8n4JNWVmc2bkzbpsNo8OB0eGgoHFjLHv2kNulS9EAqDhZB6s17Avryqu1v1Xh1pr2AwbwXsuW9Dx0CGW3Q0ICGtiUn0/PCLQqtpZjcZ0pM5MO//d/LKlXj/kFBbx5331hqFnNZ1KK1Ph4GlVgdp7JYKCTzUbLuDjW5eVVKNVKuEXnf2M15LHZjtu7VhuNbL3xRmLy8vCazewfOZKj/ftLgChDir9fN1qZDAZax8Wx2d/FU/DII5xzzTU8uGgR/UePBnzjBO0slkr3FVdGrtvNwXK80dRfuBBjbi7XuN2cd//9IclrVNvUMZnoFR9f6UkXVqORtMREDhcWsiYvr8L7hoSDDGYHCOcOd+Jk4coMW1WB6yo8bjfX9+1L/7p1uf2118g75RQAmsbG0iuMrYrAVCNleekf/+DPzZt5Mz0dazXemjQatLdaaW+xBK2bVGvNNoeDjfn5lRrwDtVgtny8FVGji80W9UECwGww0MrfxWiMieGcCRN4d8UKWtxzT1GZPU4nueUYLwgGp9dLhtNZZjnrtm2sXLSI+X/8wUV33CFBogriDAYGJiXRwWoN6liaUoq2FgtDk5NpEEXrPiRQiKjQJDaWllE2y6k0bS0WjP43iGFXXcXrMTGcsngxCevXF5VZXc5Fb1W10+EoM6dT3cWLGdGvH3vvvZd6TZowqphJF6J8GpnNnJacTN0QLq60Go30S0ykR3x80esskiRQ1EIK32K2aJFgNNKjGu0/Dce3KuKTklhzySXYgRYBq+uPuFzsLmd3UGV5tWZ7Wdfweuk8dSpZderw3tatjLvnHszVKChHC4NSdLXZ6JOYiDlMY40t4uI4LTmZlAhP7pBAUUskGI30iI9nSHIyo+vWZWSdOgxPSaFXQgLNYmMjNkUvRinSEhKiepe4kgS2KobeeiuvA63mzMG2dWtRmbV5eThDOFc+w+kscy5+06++InnlSh6Li6NOmzaMGD8+ZPWpqeIMBgYmJtI6AhMtbEYjg5KSaG+1EqmPd9Xvv1NUiEEpOlmtDElOpkVcHEkxMRiUQimFzWikaWwsPRMSGJGcTDuLBVOYA0ZqfHxYZwcFU6zBUNRd1qRtWxYMHUomYFu6tKiMS+ty512qjG1lLLAzOBx0/Oc/2dOiBf/eu5fx990XkRTi1Vkdk4khyckRzeOllKKD1cqApCTiIvChSgJFDZbg3yfhFKu1zBZDnNFIR5uNESkptLdaw9I11dlmi2hW2GAIXDB16h130ExrPjhhmuoep7NcU1cr6kBhYZlTKZNWr8Z85Ah3ezw069CBweeeG/R61GRNY2MZkJhIbJS0eOuaTEUL+sIpOn57EXTJMTEMTEqq8Cwik8FAB6uVESkpx3WtBFtbi4W2UbxeorxiA8Yqug4cSLPUVL5+5RVsa9ceV26F3Y4jyPPjN/vzTpUms08fHn30UWbt2cOl991XY3ZyC4dWcXH0jI+PupXTJoOB3gkJ9EpICFsPgASKGqiOycSAKg64mQ0GOvtbGMEOGM1iY+lktQbtfJF2iv/5UUpx/s03c962bQwdNoz4zZuLyji9Xv6y28u141x5HC4sJLOM6beJa9bgLizknVdfpXWXLgwYOzYo164N2lksdIvCNDKBmsbGMjQ5mUZhmEYrgaKGqWMy0T8xMWiDw7H+gHF6SkpQxjCax8ZGZR6nqggcqxgwdiw/Nm2KQynaPffcceWOulysL0croDw2lTE2Eb95M6eefjrqhhvYv2MHl02eHJQNe2xGI20tFgYkJjI4KYn+iYmkJSTQKi4ubDOBQq291UrHajILL85opE9iIv0TE0O6Bqlm/GUF4Nu7uW9CQki6i8wGQ9EYRiertcIDagroZrORmpAQdU35YDjWqjDGxDD45pt52eul6ZdfkrBu3XHlthUUsLcci+NKc9Tl4ojLVWqZTo89hsdiYfKSJbTv1avKe02kxMRwalISw1NS6GyzUc9sJsVkor7ZTOPYWLrFx3NGSgp9ExOpEwWbN1VWG4uFDtWwtVvfbA7p3i0RCRRKqTpKqflKqc3+7ykllPMopVb4v2aHu57VSZzBQP/EREwh/lRnMhg4xT+GkRofX67kfVajkQFJSbSqAWMSJQkcqzj90kt5NSmJXKORzo89dlLZ5XY7mWW80ZdmcxmtibqLFtFo3jy+HjCADQcOcMXDD1e6BRdrMNAzPp7ByckklxEADErR0GxmUFISfRISSKhm4yEt4+LoUk1aEsUxKBWyMcVItSgeABZordsBC/y3i1OgtU71f50dvupVLzFK0S8xMaw7wRmUorl/MdDQ5GQ6WK0kxcQQazBgUgqTUjQ2m+mfmMjwEK9ijRbHxnLibDb6XXMN09xubCtXEnvgwHHlvFqzJDe3UunIj7pcpc6gUi4XXR96iLymTbnpr79IHTqUHqeeWuHrADTwf0ptVonFeY1iYzktOZmuNltULe4sSdPYWLpV4yARapEKFOcAM/w/zwDOjVA9aoTyfrIPlYSYGNr712qcWacOI+vWZWTduqQlJlK/mM1baqpYf2ZZgLMnTeIti4XzBw7E2bDhSWULvV7+yMmp0GI8j9assNtLLWPdtQtTTg6v9+nDoaNHueKhhyr2S+DrJuxotdI3IaFK4w5KKVpbLAyLwHTOimhkNpNaw8bNgi1SgaKh1nqf/+f9wMn/ST5xSql0pdQfSqlzSzuhUmqSv2z6oUOHglnXqNbOYqn2axFqkmOtisS6dTnjuuuY/7//kbFmDfEbN55UNt/j4c+cnHLvcrYxP7/MVkhe27b893//Y8r33zPwH/+gXWpqhepvUor+iYm0C2Kyuzijkd4JCfRJSIia9QjH1DOZ6F1Dx82CKWR/NaXUD0qpNcV8nRNYTvvynJc0Z7ClP+3tpcALSqm2JV1Pa/2m1jpNa51Wv3794P0iUayB2VwtB95qssAcUOfdfDNxVittr7iCAeefj7GY1kC2281v2dllrrE46nKVuc1pozlzMDidzHztNZwOB5dPnlyhulsMBgYlJVEvRNMtG/mnczaNkg82dUwm+kiQKJeQBQqt9ela667FfH0NHFBKNQbwfz9Ywjn2+L9vAxYCPUNV3+rGajTSS5rLUamoVVGnDmOvu45HMjKIPXSITo8/Xmz5XI+HxTk55JcQLDxas7KMLqe6ixfT56qrSH76aea8+y5nXH45zdq1K3edE2NiGJyUREKIuzDNBgO9/IvFIjl2UT/I08hrukg9S7OBq/w/XwV8fWIBpVSKUirW/3M9YBCw7sRytZECesbHh3yGk6icwHUV5950E6vi4/m0VStav/sudRcvLvYx+R4Pi7Kz2ZyfX9QVpbUmw+Hg56ws7KW0OIx5eXS/6y7yWrTg3pUribNaK9SaSI6JYWBiInFhnAzR1D/YHYmsqI3NZvomJkZF+u7qIlLvNE8DZyilNgOn+2+jlEpTSr3tL9MJSFdKrQR+Ap7WWkugANpZrdV6rnptcCwHVEJKCv+YNImrt28nq0kTetxxB8YSFt0Ver1syM/nh8xMVtrt/JSVxXK7vcxxic6PPopt+3ZmXX45vy9axPj77iv3FqdJMTFhmVZdHKvRyMCkpLCmcmljsciYRCXIVqgBqsNWqMkxMQxKSpIXejWwNi+PbQUF2LOzuXHAAM6uU4cXTSbS33+fgpYtg3KNBt9/T7/LLmPT9ddz2rx5mMxm/r1wITHl+CCRFBPDgAgFiRMdKCxkeW4urhC9H5mUomdCAg2jaNe4aCNbodYQRqXoJZ+Gqo22cXEYlCI+KYkJU6bw4caNPHbddUELEgAFzZuTcf75PJOSwv4dO7jm8cfLFSQsYVqgWV4NzWZfKu8QdEXVN5k4LTlZgkQVRMerRJRLR6u1WuwpLXzijMaiGVDDLrmEDmlpvPfEEzgOHCD1lluI37Ch0udWbjd4veR26sTXd93Fh88/T//Ro+k1fHiZjzUqRd8w7tJWXlb/Bj2dbbagfBhK9Her9U9KCuti1Jooul4pokSJMTFFi7lE9XGsVWEwGLj+6afJOXyYuU89Rf2ffqLfpZdiPljshL/Seb30uOMOUm+9FbfTyfM33YQlPp4bn322XA/vFeEFmqVRStHWYmFIUlKlV/PX8a+NGJKURH1pRQSFBIpqorvNJlNhq6E4o5GW/nUDp/TowVlXXsmMTz7hy6lTiT18mL5XXlni4HaxtKbzo4/SfNYs8lu1Ytbzz7N11Spufu45Uho0KPPhHa1WGkXJOobSJPj3UxmUlFSuNNo2o5E2FgtDk5MZlJREk9hY+X8Jouj8WCGO0yIuLqLbMIqqOcViYafTiVdrLn/wQdLnz+fOZ56h2bPPcuqttzJ41CiWvvce+W3aAODxeFi2YAFbV6+mbbdu9BoxwrfhkNdL+2efpe1rr7Ht2mv5ZtgwPhs7luGXXMKAMWPKrEc9k4lTqllixjomE3VMJhweD1luNzkeD3keDyalMBsMxBkM1DOZsErXUkhJoIhyZoOhRm3yUxsdG6vYVlBAYp063PfOOzx49tnc/d//8uLHH9M9YM2Dx+Nh6kUXsWnZMpz5+cRarbTv1YvHPvuMPjfeSNOvvmL3xRfzyy238Ow//kHdxo257skny6yDSSl6VuMFmnFGI42MRhpFuiK1lASKKNfBYom6QUdRce0tFnY7HLi0pmNaGtc8/jhvPPAAr6SlMe7339ExMaA1bc85hzHLllHgchEDdM7LY316OssWLKDhpZdycMQI1p9xBlPPOYfco0f559dfY0tMLPP6PeLjw7qgTtQsEiiimM1opIUMYNcIJoOBjlYrq/PyABg9cSIb0tP5+JlnqNe0KSPGjcOUlUX8jh086XIR2Ea4u6CAbatXc/juu3Hk5THtwgvZu20bj86axSk9epR57eaxsZI4UlSJfFSNYh2tVlkzUYO0jIsr2sxHKcXN06fTddAg/n3rrcyYNg1HYiLv/etfdLBYuAG4HF9ys/esVtp068ahjAweGz+ezcuWcc8bb9B98OAyr2k1Gukq+yyIKpJAEaWSY2KiOoe/qDil1HE7qMXZbDz26aeMnDCBL156iaeuuooWnToRl5bGDJuNj5Vig81Gs169OLBrF7cMHszWlSu545VXGDh2bNnXwzcVVhLfiaqSFB4BoimFx8AqzCMX0W1pTg77A3ap01oz5913eeuhh/B6PHTs25fm7dtzZM8evF4vRw8cYNeGDaQOHcrNzz1HwxYtynWd9larpKEX5VZaCg8Zo4hCDcxmCRI1WI/4eHKys4vSiiulGHPNNfQcNoxFX33Fr7NnM/+jjwBIql+fFh06cPtLLzH8kkvKPWspJSaG9tVsKqyIXtKiCBAtLYohyckkRenKWREcuW43v2Zn4y7h/+/Arl1YbDYS69at8LlNSnFqcrKkexEVIkkBq5HGZrMEiVogISaGXgkJJR5v2KJFpYKEwZ/HSYKECCYJFFFG+pRrj4ZmM93j44M2s00BvePjZa8SEXTy0TWKNI2NDflWlCK6tIyLo77JxNq8vOMGuCujW3x8tcjjJKofeVeKEgppTdRWVqORPomJHHG5OFRYSI7HQ47bTYF/S9SyJBiNdIuPlwkQImQkUESJ5nFx0q9cy9U1mY57s3d5veR6POR6PGS73b6keG43GohRiliDgRaxsbTxb7sqRKhEJFAopS4CHsW3L3ZfrXWxU5SUUiOBfwNG4G2t9dNhq2QYxShFB5nKKE5gMhioYzAcN+bg8c+SMkpgEGEUqcHsNcD5wC8lFVBKGYFXgFFAZ2C8UqpzeKoXXqdYLJKwTZSLUSkJEiLsItKi0FqvB8paPNQX2KK13uYv+wlwDrAu5BUMI4vBQFtpTQgholg0T49tCuwOuJ3hv69YSqlJSql0pVT6oUOHQl65YAnW/sBCCBEqIWtRKKV+gGL3GXlIa/11sK+ntX4TeBN8K7ODff5QqGMySeI/IUTUC1mg0FqfXsVT7AGaB9xu5r+vRjAoRXdJ/yyEqAaiuetpKdBOKdVaKWUGxgGzI1ynoOlitcriOiFEtRCRQKGUOk8plQEMAL5VSs3z399EKTUHQGvtBm4B5gHrgU+11msjUd9ga2g200oGsIUQ1USkZj19BXxVzP17gdEBt+cAc8JYtZCLNRhIjY+PdDWEEKLcornrqcYxKEWv+HjMsuOYEKIakXesMOoZH089sznS1RBCiAqRQBEm3Ww2mQorhKiWJFCEQXurVQavhRDVlszPDLEuNhttJEgIIaoxCRQhYlCKnvHx0t0khKj2JFCEQKzBQO+EBNlIRghRI0igCLJ6JhM94+MlbbgQosaQQBEkCmhntdLeYikrfboQQlQrEiiCoK7JRDebTXI3CSFqJHlnqwKLwUBHq5VmcXGRrooQQoSMBIpKiDcaOcVioWlsrGw6JISo8SRQVEBiTAztLBYam80yDiGEqDUkUJSDzWiks9VKQwkQQohaSAJFGZrExtLDZiNGMr4KIWopCRQlMChFV5uNljJQLYSo5SRQlKC7zUZzCRJCCCHZY4vTKi5OgoQQQvhFas/si5RSa5VSXqVUWinldiilViulViil0sNRt5SYGLrYbOG4lBBCVAuR6npaA5wPvFGOssO01odDXB8ATAYDaQkJsjZCCCECRCRQaK3XA1E31dQmifyEEOIk0T5GoYHvlVJ/KaUmRboyQghRG4WsRaGU+gFoVMyhh7TWX5fzNIO11nuUUg2A+UqpDVrrX0q43iRgEkCLFi0qVWchhBAnC1mg0FqfHoRz7PF/P6iU+groCxQbKLTWbwJvAqSlpemqXlsIIYRP1HY9KaVsSqmEYz8DZ+IbBBdCCBFGkZoee55SKgMYAHyrlJrnv7+JUmqOv1hD4Fel1EpgCfCt1vq7SNRXCCFqs0jNevoK+KqY+/cCo/0/bwN6hLlqQgghThC1XU9CCCGigwQKIYQQpZJAIYQQolQSKIQQQpRKAoUQQohSSaAQQghRKgkUQgghSiWBQgghRKkkUAghhCiVBAohhBClkkAhhBCiVBIohBBClEoChRBCiFJJoBBCCFEqCRRCCCFKJYFCCCFEqSRQCCGEKJUECiGEEKWSQCGEEKJUEQkUSqlnlVIblFKrlFJfKaWSSyg3Uim1USm1RSn1QJirKYQQgsi1KOYDXbXW3YFNwOQTCyiljMArwCigMzBeKdU5rLUUQggRmUChtf5ea+323/wDaFZMsb7AFq31Nq11IfAJcE646iiEEMInJtIVACYCs4q5vymwO+B2BtCvpJMopSYBk/w37UqpjZWsTz3gcCUfG0pSr4qRelWM1KtiamK9WpZ0IGSBQin1A9ComEMPaa2/9pd5CHAD/6nq9bTWbwJvVvU8Sql0rXVaVc8TbFKvipF6VYzUq2JqW71CFii01qeXdlwpNQEYC4zQWutiiuwBmgfcbua/TwghRBhFatbTSOA+4GytdX4JxZYC7ZRSrZVSZmAcMDtcdRRCCOETqVlPLwMJwHyl1Aql1OsASqkmSqk5AP7B7luAecB64FOt9dow1K3K3VchIvWqGKlXxUi9KqZW1UsV3+sjhBBC+MjKbCGEEKWSQCGEEKJUtTJQKKUuUkqtVUp5lVIlTiUrKYWIf4D9T//9s/yD7cGoVx2l1Hyl1Gb/95Riygzzj+sc+3Iopc71H3tfKbU94FhquOrlL+cJuPbsgPsj+XylKqV+9/+9VymlLgk4FtTnq6yUM0qpWP/vv8X/fLQKODbZf/9GpdRZValHJep1l1Jqnf/5WaCUahlwrNi/aZjqNUEpdSjg+tcGHLvK/3ffrJS6Ksz1ej6gTpuUUlkBx0LyfCml3lVKHVRKrSnhuFJKveiv8yqlVK+AY1V/rrTWte4L6AR0ABYCaSWUMQJbgTaAGVgJdPYf+xQY5//5deDGINXr/4AH/D8/ADxTRvk6wFHA6r/9PnBhCJ6vctULsJdwf8SeL6A90M7/cxNgH5Ac7OertNdLQJmbgNf9P48DZvl/7uwvHwu09p/HGMZ6DQt4Dd14rF6l/U3DVK8JwMvFPLYOsM3/PcX/c0q46nVC+VuBd8PwfA0BegFrSjg+GpgLKKA/8Gcwn6ta2aLQWq/XWpe1crvYFCJKKQUMBz73l5sBnBukqp3jP195z3shMFeXPMU4WCparyKRfr601pu01pv9P+8FDgL1g3T9QOVJORNY38+BEf7n5xzgE621U2u9HdjiP19Y6qW1/ingNVRSSp1gq0qKnrOA+Vrro1rrTHy540ZGqF7jgZlBunaJtNa/4PtQWJJzgA+0zx9AslKqMUF6rmploCin4lKINAXqAln671xVx+4PhoZa633+n/cDDcsoP46TX6T/9Dc9n1dKxYa5XnFKqXSl1B/HusOIoudLKdUX36fErQF3B+v5Kun1UmwZ//ORje/5Kc9jQ1mvQNfg+2R6THF/03DW6wL/3+dzpdSxBbhR8Xz5u+haAz8G3B2q56ssJdU7KM9VNOR6CglVjhQikVBavQJvaK21UqrEucv+Twvd8K0zOWYyvjdMM7751PcD08JYr5Za6z1KqTbAj0qp1fjeDCstyM/Xh8BVWmuv/+5KP181kVLqciANOC3g7pP+plrrrcWfIei+AWZqrZ1KqevxtcaGh+na5TEO+Fxr7Qm4L5LPV8jU2EChy0ghUg4lpRA5gq9ZF+P/VFih1CKl1UspdUAp1Vhrvc//xnawlFNdDHyltXYFnPvYp2unUuo94J5w1ktrvcf/fZtSaiHQE/iCCD9fSqlE4Ft8HxL+CDh3pZ+vYpQn5cyxMhlKqRggCd/rKZTpasp1bqXU6fiC72laa+ex+0v4mwbjja/MemmtjwTcfBvfmNSxxw494bELg1CnctUrwDjg5sA7Qvh8laWkegfluZKup5IVm0JE+0aIfsI3PgBwFRCsFsps//nKc96T+kb9b5bHxgXOBYqdIRGKeimlUo513Sil6gGDgHWRfr78f7uv8PXffn7CsWA+X+VJORNY3wuBH/3Pz2xgnPLNimoNtAOWVKEuFaqXUqon8Aa+lDoHA+4v9m8axno1Drh5Nr4MDeBrRZ/pr18KcCbHt6xDWi9/3TriGxz+PeC+UD5fZZkNXOmf/dQfyPZ/EArOcxWKEfpo/wLOw9dX5wQOAPP89zcB5gSUG41vY6Wt+D6NHru/Db5/5C3AZ0BskOpVF1gAbAZ+AOr4708D3g4o1wrfJwXDCY//EViN7w3vIyA+XPUCBvqvvdL//ZpoeL6AywEXsCLgKzUUz1dxrxd8XVln+3+O8//+W/zPR5uAxz7kf9xGYFSQX+9l1esH///Bsedndll/0zDV6ylgrf/6PwEdAx470f88bgGuDme9/LcfBZ4+4XEhe77wfSjc538tZ+AbS7oBuMF/XOHb6G2r/9ppAY+t8nMlKTyEEEKUSrqehBBClEoChRBCiFJJoBBCCFEqCRRCCCFKJYFCCCFEqSRQCCGEKJUECiGEEKWSQCFEGCilflJKneH/+Qml1EuRrpMQ5VVjcz0JEWWmAtOUUg3w5f85O8L1EaLcZGW2EGGilPoZiAeGaq1zI10fIcpLup6ECAOlVDegMVAoQUJUNxIohAgxfxbU/+DbhcyulArWbmxChIUECiFCSCllBb4E7tZarwcexzdeIUS1IWMUQgghSiUtCiGEEKWSQCGEEKJUEiiEEEKUSgKFEEKIUkmgEEIIUSoJFEIIIUolgUIIIUSp/h/fgVkZCNXRKwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAShElEQVR4nO3de5CddX3H8fdHIlJF5ZKVYgCDI1Cp6MikFKv1QrRVqAQr4+CoRCfTjHettoraGe1lptCpWmmtNhU1WkGQosRrq1wGb0SDUC7BS8AgwZCsVaiXeqF++8d54qxhN3t2z+6ePb+8XzM7e57n+Z3n+f6yZz/57e95znNSVUiS2nKfYRcgSZp7hrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3LUoJLkpyZOHXcd8SnJMkuuS/DDJK+don7cneexc7EttMdw175JsTfLU3da9MMkXdi1X1W9X1ZUz3c+IeR1wRVU9sKrOHXRnSQ4EDgVuHrgyNcdwl4AkSxbgMA8Dbprpk/ZQ23HAlqr66UBVqUmGuxaFiaPyJK9Pckc3ffGNJCuTfBA4Avh4kh8leV3X9pFJrkxyVze1c+qEfR6f5NpuPx9JcmGSv9ntmK9Pcj3w4yRLkpyV5JbuOZuTPGu39n+e5PokP05yXpJDkny6a/+5bjQ9Wf8uB54C/FNX/9HT1H6v2ibZ7aOBG7v2909yfpJLkuw/+5+EWmG4a1FJcgzwcuB3quqBwB8CW6vqBcB3gGdW1f5V9XdJ7gt8HPhP4CHAK4APdXPb+wIfBd4PHARcADzrXgeE5wKnAAdU1T3ALcDvAw8G/hL4tySHTmj/bOBpwNHAM4FPA28Exuj9Pk06l15VJwGfB15eVfsD356q9j3UtrvjgBuSHAl8EfgG8Oyq+tFkNWjvshB/ikoAH0syMaD2Bb42Sbv/A+4HHJtkvKq27mGfJwL7A2dX1S+By5N8gl4oXk7v9X1u9W59ekmSr0yyj3Or6vZdC1X1kQnbLkzyBuAE4NJu3T9W1Q6AJJ8HdlbVtd3yR4GVe6i339rfMlltk3g0UMAVwKuq6tI9tNVexpG7FsppVXXAri/gpZM1qqotwKvpBdzOJB9O8tAp9vlQ4PYuHHe5DVjWbbujfv2e1pMF5a+tS3Jmd0XLXUnuAh4FLJ3QZMeEx/87yXK/UyJ7qn1P9e6qM11tzwLeNVWwJ/F3fC/lD16LTlWdX1VPoHcCsoBzdm3arel3gcN3C7AjgDuA7cCyLgR3OXyyw+16kORhwL/SmxY6uPtP6EYgkzxvUHuq/V61TeLI7vtTgdcmWTFxY5JrkrybXn+0FzLctah08+UnJbkf8FN6o+Fdo9sdwMMnNN8I/AR4XZL7dtfJPxP4MPBlelM8L+9OlK6iN72yJw+gF6jjXS0vojc6ng97qr0fjwaur6obgLXAR3edG0iylN48/huras1cF67RYLhrsbkfcDbwPeBOeiH1hm7b3wJ/0U2Z/FlV/ZxeID6ja//PwJlV9fVu2x8Da4C7gOcDnwB+NtWBq2oz8FZ6/zHsoHfC8otz3cHuWFPW3ucujgOu7/b1MWAdvfMa+9EL/vOr6vtzXbdGR/yYPe0tkmwE3l1V7xt2LfMpyauBbVV18bBr0fA4clezkjwpyW920zKr6Y1oPzPsuhbAccB1wy5Cw+WlkGrZMcBF9ObSbwVOr6rtwy1p/jnPLnBaRpKa5LSMJDVoUUzLLF26tJYvXz7sMiRppFxzzTXfq6qxybYtinBfvnw5mzZtGnYZkjRSktw21TanZSSpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGL4h2qo2r5WZ/sq93Ws0+Z50ok6dc5cpekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQtOGe5L1Jdia5ccK6g5J8Nsm3uu8HduuT5NwkW5Jcn+T4+SxekjS5fkbu7weevtu6s4DLquoo4LJuGeAZwFHd11rgXXNTpiRpJqYN96q6Cvj+bqtXAeu7x+uB0yas/0D1XA0ckOTQOapVktSn2c65H1JV27vHdwKHdI+XAbdPaLetW3cvSdYm2ZRk0/j4+CzLkCRNZuATqlVVQM3ieeuqakVVrRgbGxu0DEnSBLMN9x27plu67zu79XcAh09od1i3TpK0gGYb7huA1d3j1cClE9af2V01cyJw94TpG0nSApn2A7KTXAA8GViaZBvwZuBs4KIka4DbgOd0zT8FnAxsAX4CvGgeapYkTWPacK+q506xaeUkbQt42aBFSZIG4ztUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg5YMuwDt3Zaf9cm+2m09+5R5rkRqiyN3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KCBwj3Jnya5KcmNSS5Isl+SI5NsTLIlyYVJ9p2rYiVJ/Zl1uCdZBrwSWFFVjwL2Ac4AzgHeXlWPAH4ArJmLQiVJ/Rt0WmYJ8BtJlgD3B7YDJwEXd9vXA6cNeAxJ0gzNOtyr6g7g74Hv0Av1u4FrgLuq6p6u2TZg2WTPT7I2yaYkm8bHx2dbhiRpEoNMyxwIrAKOBB4KPAB4er/Pr6p1VbWiqlaMjY3NtgxJ0iQGmZZ5KvDtqhqvql8AlwCPBw7opmkADgPuGLBGSdIMDRLu3wFOTHL/JAFWApuBK4DTuzargUsHK1GSNFODzLlvpHfi9GvADd2+1gGvB16TZAtwMHDeHNQpSZqBgW75W1VvBt682+pbgRMG2a8kaTC+Q1WSGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0a6AOypaksP+uTwy5B2qs5cpekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1aKBwT3JAkouTfD3JzUkel+SgJJ9N8q3u+4FzVawkqT+DjtzfAXymqn4LeAxwM3AWcFlVHQVc1i1LkhbQrMM9yYOBJwLnAVTVz6vqLmAVsL5rth44bbASJUkzNcjI/UhgHHhfkmuTvCfJA4BDqmp71+ZO4JDJnpxkbZJNSTaNj48PUIYkaXeDhPsS4HjgXVX1WODH7DYFU1UF1GRPrqp1VbWiqlaMjY0NUIYkaXeDhPs2YFtVbeyWL6YX9juSHArQfd85WImSpJmadbhX1Z3A7UmO6VatBDYDG4DV3brVwKUDVShJmrFB7+f+CuBDSfYFbgVeRO8/jIuSrAFuA54z4DEkSTM0ULhX1XXAikk2rRxkv5KkwfhJTJrRpyZtPfuUeaxkav3WOKz6pMXG2w9IUoMMd0lqkNMyDfNDqqW9lyN3SWqQ4S5JDTLcJalBhrskNchwl6QGebXMCPIqGEnTceQuSQ0y3CWpQYa7JDXIOfdFxLn0wXmDManHkbskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkJdCaka8XFMaDY7cJalBhrskNchwl6QGjfycu28313yayTkGX2NaTBy5S1KDRn7kPgq8wkTSQnPkLkkNMtwlqUEDh3uSfZJcm+QT3fKRSTYm2ZLkwiT7Dl6mJGkm5mLk/irg5gnL5wBvr6pHAD8A1szBMSRJMzBQuCc5DDgFeE+3HOAk4OKuyXrgtEGOIUmauUFH7v8AvA74Zbd8MHBXVd3TLW8Dlk32xCRrk2xKsml8fHzAMiRJE8063JP8EbCzqq6ZzfOral1VraiqFWNjY7MtQ5I0iUGuc388cGqSk4H9gAcB7wAOSLKkG70fBtwxeJmSpJmY9ci9qt5QVYdV1XLgDODyqnoecAVwetdsNXDpwFVKkmZkPq5zfz3wmiRb6M3BnzcPx5Ak7cGc3H6gqq4Eruwe3wqcMBf7lSTNju9QlaQGeeOwSXijL0mjzpG7JDXIcJekBhnuktQg59y1V/K8ilrnyF2SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBe81dIb0LoKS9iSN3SWqQ4S5JDdprpmWkxaLfKcKtZ58yz5WoZY7cJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoNmHe5JDk9yRZLNSW5K8qpu/UFJPpvkW933A+euXElSPwYZud8DvLaqjgVOBF6W5FjgLOCyqjoKuKxbliQtoFmHe1Vtr6qvdY9/CNwMLANWAeu7ZuuB0wasUZI0Q3PyDtUky4HHAhuBQ6pqe7fpTuCQKZ6zFlgLcMQRR8xFGVJTfCerBjHwCdUk+wP/Dry6qv5n4raqKqAme15VrauqFVW1YmxsbNAyJEkTDDRyT3JfesH+oaq6pFu9I8mhVbU9yaHAzkGLlEaBt5XWYjLI1TIBzgNurqq3Tdi0AVjdPV4NXDr78iRJszHIyP3xwAuAG5Jc1617I3A2cFGSNcBtwHMGqlDSHs313Lxz/W2YdbhX1ReATLF55Wz3K0kanO9QlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJatCcfMyepL3PTD6cxNsDLzxH7pLUIEfukuadHwCy8By5S1KDHLlLewk/wHvv4shdkhrkyF3SyHEOf3qGu6RFw6mjueO0jCQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDfI6d0nNmuvr5kfpTVHzMnJP8vQk30iyJclZ83EMSdLU5nzknmQf4J3A04BtwFeTbKiqzXN9LElarIZ9i4T5GLmfAGypqlur6ufAh4FV83AcSdIU5mPOfRlw+4TlbcDv7t4oyVpgbbf4oyTfmOXxlgLfm+VzF4NRrx9Gvw+jXj+Mfh9Gov6cs8fNs+rDNPuczsOm2jC0E6pVtQ5YN+h+kmyqqhVzUNJQjHr9MPp9GPX6YfT7MOr1w+Lrw3xMy9wBHD5h+bBunSRpgcxHuH8VOCrJkUn2Bc4ANszDcSRJU5jzaZmquifJy4H/APYB3ltVN831cSYYeGpnyEa9fhj9Pox6/TD6fRj1+mGR9SFVNewaJElzzNsPSFKDDHdJatDIhPt0tzRIcr8kF3bbNyZZPoQyp9RH/a9JsjnJ9UkuSzLl9avD0u9tJZI8O0klWTSXhUF/9Sd5TvdzuCnJ+Qtd43T6eB0dkeSKJNd2r6WTh1HnZJK8N8nOJDdOsT1Jzu36dn2S4xe6xun00YfndbXfkORLSR6z0DX+SlUt+i96J2ZvAR4O7Av8F3Dsbm1eCry7e3wGcOGw655h/U8B7t89fsliqr/fPnTtHghcBVwNrBh23TP8GRwFXAsc2C0/ZNh1z6IP64CXdI+PBbYOu+4JtT0ROB64cYrtJwOfBgKcCGwcds2z6MPvTXj9PGOYfRiVkXs/tzRYBazvHl8MrEySBaxxT6atv6quqKqfdItX03t/wGLS720l/ho4B/jpQhbXh37q/xPgnVX1A4Cq2rnANU6nnz4U8KDu8YOB7y5gfXtUVVcB399Dk1XAB6rnauCAJIcuTHX9ma4PVfWlXa8fhvx7PCrhPtktDZZN1aaq7gHuBg5ekOqm10/9E62hN4JZTKbtQ/dn9OFVNbf3WZ0b/fwMjgaOTvLFJFcnefqCVdeffvrwFuD5SbYBnwJesTClzYmZ/p4sdkP9PfZ+7otMkucDK4AnDbuWmUhyH+BtwAuHXMogltCbmnkyvRHXVUmOq6q7hlnUDD0XeH9VvTXJ44APJnlUVf1y2IXtTZI8hV64P2FYNYzKyL2fWxr8qk2SJfT+JP3vBaluen3dkiHJU4E3AadW1c8WqLZ+TdeHBwKPAq5MspXenOmGRXRStZ+fwTZgQ1X9oqq+DXyTXtgvFv30YQ1wEUBVfRnYj94NrUZBE7cuSfJo4D3AqqoaWgaNSrj3c0uDDcDq7vHpwOXVndVYBKatP8ljgX+hF+yLba4XpulDVd1dVUuranlVLac333hqVW0aTrn30s9r6GP0Ru0kWUpvmubWBaxxOv304TvASoAkj6QX7uMLWuXsbQDO7K6aORG4u6q2D7uomUhyBHAJ8IKq+uZQixn22ecZnKU+md5I6hbgTd26v6IXINB7EX8E2AJ8BXj4sGueYf2fA3YA13VfG4Zd80z7sFvbK1lEV8v0+TMIvamlzcANwBnDrnkWfTgW+CK9K2muA/5g2DVPqP0CYDvwC3p/Ja0BXgy8eMK//zu7vt2w2F4/ffbhPcAPJvwebxpWrd5+QJIaNCrTMpKkGTDcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoP+H5LbJ+rSoxSaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.41312072 0.24882542\n"
     ]
    }
   ],
   "source": [
    "u_pred, k_pred = model.predict(x_test, samples, [process_u, process_k])\n",
    "plot1d(x_u_train, u_train, x_test, u_test, u_pred[..., 0], ylabel=\"$u$\", ylim=[-2.0, 2.0], title=\"$HalfN(0, 1^2)$\")\n",
    "plt.hist(k_pred, bins=30)\n",
    "plt.title(\"Histogram for $k_r$\")\n",
    "plt.show()\n",
    "print(np.mean(k_pred), np.std(k_pred))\n",
    "\n",
    "# sio.savemat(\n",
    "#     \"HalfNormal.mat\",\n",
    "#     {\n",
    "#         \"x_u_train\": x_u_train, \"u_train\": u_train,\n",
    "#         \"x_f_train\": x_f_train, \"f_train\": f_train,\n",
    "#         \"x_test\": x_test, \"u_test\": u_test, \"f_test\": f_test,\n",
    "#         \"u_pred\": u_pred, \"f_pred\": f_pred, \"k_pred\": k_pred,\n",
    "#         \"noise\": noise,\n",
    "#     }\n",
    "# )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af068643",
   "metadata": {},
   "source": [
    "#### Log Normal prior for $k_r$\n",
    "\n",
    "$X\\sim LogN(0, 1^2)$ is equivalent to $\\log X\\sim N(0, 1^2)$. Choosing this prior for $k_r$ indicates that we at least known $k_r> 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ff1a8c52",
   "metadata": {},
   "outputs": [],
   "source": [
    "class LogNormal(neuq.variables._Samplable):\n",
    "    def __init__(self, mean, sigma=1.0, shape=[], initializer=None):\n",
    "        super().__init__()\n",
    "\n",
    "        self._num_tensors = 1\n",
    "        # initial value should be set to a positive value\n",
    "        self._initial_values = [tf.constant(1.0, tf.float32)]\n",
    "        if shape == []:\n",
    "            # make sure constant has at least 1 dimension\n",
    "            self._initial_values = [self.initial_values[0][None, ...]]\n",
    "        self.dist = tfp.distributions.LogNormal(loc=mean, scale=sigma)\n",
    "\n",
    "    def log_prob(self, samples):\n",
    "        # Note: here, because a constant is considered, `samples` is a list of only\n",
    "        # one element.\n",
    "        return self.dist.log_prob(samples[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9cfb04a0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Supporting backend tensorflow.compat.v1\n",
      "\n",
      "Compiling a MCMC method\n",
      "\n",
      "sampling from posterior distribution ...\n",
      "\n",
      "Finished sampling from posterior distribution ...\n",
      "\n",
      "Acceptance rate:  0.669\n"
     ]
    }
   ],
   "source": [
    "prior_k = LogNormal(mean=0, sigma=1)\n",
    "\n",
    "process_k = neuq.process.Process(\n",
    "    surrogate=neuq.surrogates.Identity(),\n",
    "    prior=prior_k,\n",
    ")\n",
    "samples, model = build_model(process_u, process_k, burnin_steps=2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "5aafc07d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEaCAYAAAAPGBBTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABWtUlEQVR4nO3dd3hTZfvA8e+TpCPpZm8Zlj0KlCEgL4gKiCK+DkB9FVFxgLgHiqLgwI0b+SluEdyo4EIQBAUqe8kSoezR3SbNeH5/JNRQupvV9v5cV68m55yc8zRtc59n3Y/SWiOEEEIUxxDsAgghhAhtEiiEEEKUSAKFEEKIEkmgEEIIUSIJFEIIIUokgUIIIUSJJFAIIYQokQQKIaoJpVRPpdTvSqmlSqk5SqmwYJdJVA8SKISoPvYB52it+wN7gIuDWxxRXUigECJAlFJPKaXu8Nf5tdYHtdZ5nqf5gMvr2quUUh38dW1RvUmgEDWKUqqOUkorpRr78JwJnnP+Xmj7TKXUi57HdYFrgDe99tdSSn2plMpRSv2jlLqyjNeboJRKUUrZlFLvFrH/DOB84Buvzc8BU8v9wwmBBApR8yQBR7XW+318zkNAe6VUA6/tXYF1nsdjgAVed/wAr+G+868PXAW8Uca7/gPA48DswjuUUrHAB8AYrbXda9d8YGCh8glRJhIoRE2TBKwvvFG53e+5s09XSs1TSsV57TcopR5RSh1XSh1QSo1WSuUrpRI850wBfsLTL6CUMgKdgLWeUwwFfvU6XxRwKfCw1jpba/0b7g/z/5X2A2itv9BafwUcL/QzmIBPgMe01n8Veo0V+BMYXNr5hShMAoWoaZL49y7f2zTcH+a9gQZABPCI1/5HgXOBzkB74A7gsNY6jX9rDl8BIzzHt8X9/7XV87wT4P3h3RpwaK23e21bD1SmH2E00At4WCm1RCk1stD+rUCXSpxf1FCmYBdAiABLAhZ6b1BK1QduA9pprQ96tn0G3OB5XBd3YOh2sslKKfU90MPrnF8DvwAzlVIxnm2bvZp/4oEsr8tGA5mFypYBxFT0B9Naf4C72ak4WUDDip5f1FxSoxA1hlIqEmjDv81BJ50NbNRaH/DaVgc46Hk8CNigtd7ptb8WsFEpFQG0A9Z5ahercNdMvPsnANI4NQhkA7GFyhHLqcHE12KAdD+eX1RTEihETdIRsHNqExBAXdx3894uBn7zPK6D1wesp/9hKLDBc85cYLdn91e4m5+6cmpA2oC7uemk7YBJKZXota0LsLnsP065taOI/hkhSiOBQtQkScAWIEwpFen5igBWA2cppVoppaKVUlNxj0Q6OapoK9BXKXWmZ1TRy0ArYCPugLBB/7tU5HzgAk6vUSwA/nPyidY6B/gCmKqUilJK9cUdnD4AUEq9W9TQV88+k6d2ZASMnp+jxGZkz/HdcXe4C1EuEihETZKE+8Myz+tri9Y6BXgCdw0iFfed9zla61wArfUiYA7uGkIK7gCQC2yjUOe41noP7lnR8Zx69/4+cIFSyuy17VbADBzxnP8WrfXJGkVTYHkxP8dkT9kfAK72PJ5cys9+EbCkUPOaEGWiZM1sIcpHKXUzMExrfVE5X/ckcERrPaOU48JxB5nOheZCVJhSaiVwvdZ6ky/OJ2oWCRRClEIp1Rt3x/Y+3B3bHwEXaa1XBrVgQgSIDI8VonRdge+AMNyd0GMkSIiaJGg1CqVUU9zttvUBDczSWr9U6BgFvIS7czAX9z/omkCXVQgharJg1igcwN1a6zWeCUp/KqV+0lpv8TpmKJDo+eoFvOH5LoQQIkCCNurJkxJ5jedxFu4hiIUzel4MvK/d/gDilVIys1QIIQIoJPoolFLNcbcDF273bYy7A/GkVM+2g4WOQyk1DhgHEBUV1b1t27Z+KasQQlRHf/755zGtdd2i9gU9UCilooHPgTu01oVz35SZ1noWMAsgOTlZp6Sk+KiEQghR/Sml/iluX1An3HnW9P0c+Ehr/UURh+zHPfHopCaebUIIIQIkaIHCM6LpbWCr1vqFYg6bD1zjWSugN5BxMrunEEKIwAhm01Nf3Iu0bFRKrfNsexBoBqC1nok7P84FwE7cw2OvC3wxhRCiZgtaoPCs6KVKOUYD4wNTIiGEP9jtdlJTU7FarcEuigAiIyNp0qQJYWFhZX5N0DuzhRDVW2pqKjExMTRv3hx3i7MIFq01x48fJzU1lRYtWpT5dZI9VgjhV1arldq1a0uQCAFKKWrXrl3u2p0ECiGE30mQCB0V+V1IoBBCCFEiCRRCiGrt+PHjJCUlkZSURIMGDWjcuHHB8/z8/BJfm5KSwsSJE0u9Rp8+fXxV3FMMGDCA0iYPz5gxg9zcXL9c/yTpzBZChBSn08nChQtZu3YtXbt2ZejQoRiNxgqfr3bt2qxbtw6ARx99lOjoaO65556C/Q6HA5Op6I/C5ORkkpOTS73GihUrKly+ypoxYwZXX301FovFb9eQGoUQImQ4nU4GDx7M6NGjmTJlCqNHj2bw4ME4nU6fXmfMmDHcfPPN9OrVi/vuu49Vq1Zx1lln0bVrV/r06cNff/0FwJIlS7jwwgsBd5AZO3YsAwYMoGXLlrz88ssF54uOji44fsCAAVx22WW0bduWq666ipNLOSxYsIC2bdvSvXt3Jk6cWHBeb3l5eYwaNYp27dpxySWXkJeXV7DvlltuITk5mQ4dOjBlyhQAXn75ZQ4cOMDAgQMZOHBgscdVltQohBAhY+HChaxcuZLs7GwAsrOzWblyJQsXLizyg7UyUlNTWbFiBUajkczMTJYtW4bJZOLnn3/mwQcf5PPPPz/tNdu2bWPx4sVkZWXRpk0bbrnlltPmI6xdu5bNmzfTqFEj+vbty/Lly0lOTuamm25i6dKltGjRgtGjRxdZpjfeeAOLxcLWrVvZsGED3bp1K9j3xBNPUKtWLZxOJ4MGDWLDhg1MnDiRF154gcWLF1OnTp1ij+vcuXOl3iupUQghQsbatWvJyck5ZVtOTk5B05EvXX755QVNWhkZGVx++eV07NiRO++8k82bNxf5mmHDhhEREUGdOnWoV68ehw8fPu2Ynj170qRJEwwGA0lJSezZs4dt27bRsmXLgrkLxQWKpUuXcvXVVwPQuXPnUz7g582bR7du3ejatSubN29my5YtRZ6jrMeVhwQKIUTI6Nq1K1FRUadsi4qKIikpyefX8r7Oww8/zMCBA9m0aRPffPNNsfMMIiIiCh4bjUYcDkeFjimvv//+m+eee45FixaxYcMGhg0bVmQZy3pceUmgEEKEjKFDh9KrVy+io6NRShEdHU2vXr0YOnSoX6+bkZFB48buddPeffddn5+/TZs27N69mz179gAwd+7cIo/r378/H3/8MQCbNm1iw4YNAGRmZhIVFUVcXByHDx9m4cKFBa+JiYkhKyur1OMqQ/oohBAhw2g08sMPP7Bw4ULWrVtHUlJSpUc9lcV9993Htddey+OPP86wYcN8fn6z2czrr7/OkCFDiIqKokePHkUed8stt3DdddfRrl072rVrR/fu3QHo0qULXbt2pW3btjRt2pS+ffsWvGbcuHEMGTKERo0asXjx4mKPqwx1ske+OpGFi4QIHVu3bqVdu3bBLkbQZWdnEx0djdaa8ePHk5iYyJ133hmUshT1O1FK/am1LnIssDQ9CSFEAPzf//0fSUlJdOjQgYyMDG666aZgF6nMpOlJCCEC4M477wxaDaKypEbhxerjST1CCFEdSKDwYnW5yPLBUDYhhKhOJFAUcsRuD3YRhBAipAQ1UCilZiuljiilNhWzf4BSKkMptc7z9Yi/y3SklGySQghR0wS7RvEuMKSUY5ZprZM8X1P9XaATDgcOl8vflxFCBJDRaCxILZ6UlMT06dN9du5169axYMECn50vFAV11JPWeqlSqnkwy1CYS2uO2e008JqGL4So2sxms1/yRYE7UKSkpHDBBRf45fyhINg1irI4Sym1Xim1UCnVIRAXlH4KIaq/jIwM2rRpU5BSfPTo0fzf//0fUHyq7tWrV9OnTx+6dOlCz549ycjI4JFHHmHu3LkkJSUVm5qjqgv1eRRrgDO01tlKqQuAr4DEog5USo0DxgE0a9asUhc9KoFCCL+44447fH5nn5SUxIwZM0o8Ji8v75TEgpMmTWLkyJG8+uqrjBkzhttvv520tDRuvPFGoOhU3W3btmXkyJHMnTuXHj16kJmZicViYerUqaSkpPDqq6/69OcKJSEdKLTWmV6PFyilXldK1dFaHyvi2FnALHCn8KjMdXOdTrIdDqKLWfVKCFG1FNf0dN555/Hpp58yfvx41q9fX7B93rx5zJo1C4fDwcGDB9myZQtKKRo2bFiQpyk2NjZQxQ+6kP4kVEo1AA5rrbVSqifuprLjgbj2EbtdAoUQPlbanX+guVwutm7disViIS0tjSZNmhSk6l69ejUJCQmMGTPGJ6m6q7JgD4+dA/wOtFFKpSqlrldK3ayUutlzyGXAJqXUeuBlYJQOUBZDGSYrRPX34osv0q5dOz7++GOuu+467HZ7sam627Rpw8GDB1m9ejUAWVlZOByOU9J8V1fBHvVU9DJP/+5/FQhKw99xhwO7y0WYoSr09wshSlK4j2LIkCFcd911vPXWW6xatYqYmBj69+/P448/zmOPPVZkqu7w8HDmzp3LbbfdRl5eHmazmZ9//pmBAwcyffp0kpKSCvo+qhtJM+4l3W5nWUZGwfNuMTE0lmGyQlSKpBkPPZJm3IcO2mzBLoIQQgSdBIoSHLHbcVbDGpcQQpSHBIoSOLWWTm0hRI0ngaIUByVQCCFqOAkUpTicn49Lmp+EEDWYBIpSOLSWlB5CiBpNAkUZyOgnIao2pRRXX311wXOHw0HdunW58MILg1iqqkMCRRkckuYnIaq0qKgoNm3aRF5eHgA//fQTjRs3DnKpqg4JFGVgl9FPQlR5F1xwAd999x0Ac+bMYfTofxND5OTkMHbsWHr27EnXrl35+uuvAdizZw9nn3023bp1o1u3bqxYsQKAJUuWMGDAAC677DLatm3LVVddRXWcvHySZL0ro302myxmJIQvDBhw+rYrroBbb4XcXChqAaAxY9xfx47BZZedum/JkjJddtSoUUydOpULL7yQDRs2MHbsWJYtWwa404qfc845zJ49m/T0dHr27Mm5555LvXr1+Omnn4iMjGTHjh2MHj2ak1kf1q5dy+bNm2nUqBF9+/Zl+fLl9OvXr8xvQ1UigaKMDufnY3O5iJDcT0JUSZ07d2bPnj3MmTPntNXofvzxR+bPn89zzz0HgNVqZe/evTRq1IgJEyawbt06jEYj27dvL3hNz549adKkCeBeE2PPnj0SKGo6Dey32WhpNge7KEJUbSXVACyWkvfXqVPmGkRRhg8fzj333MOSJUs4fvzfFQu01nz++ee0adPmlOMfffRR6tevz/r163G5XERGRhbsi/BqYTAajTgcjgqXK9TJ7XE57JPRT0JUaWPHjmXKlCl06tTplO2DBw/mlVdeKehnWLt2LeBeLrVhw4YYDAY++OADnE5nwMscCiRQlEOmw0FmNb5rEKK6a9KkCRMnTjxt+8MPP4zdbqdz58506NCBhx9+GIBbb72V9957jy5durBt2zaioqICXeSQIGnGvRROM16UlmYzHWroH4sQFSFpxkOPpBn3s1SbTeZUCCFqFAkU5ZTvckmiQCFEjRLsNbNnK6WOKKU2FbNfKaVeVkrtVEptUEp1C3QZi/K3Z3anEELUBMGuUbwLDClh/1Ag0fM1DngjAGUqVZrDQZokChRC1BBBDRRa66XAiRIOuRh4X7v9AcQrpRoGpnQl+9tqDXYRhBDiFP4anBTsGkVpGgP7vJ6nerYF3QGbDWsNHVMthAg9DpcLl5/OXW1mZiulxuFunqJZs2Z+v54G/rHZaGOx+P1aQlQn3xw75tPzXVSnTon79+3bxzXXXMPhw4dRSjFu3Dhuv/32Mp173bp1HDhw4LSUHyc1b96clJQU6pRShkDI15oIpfxy7lCvUewHmno9b+LZdhqt9SytdbLWOrlu3boBKdweq1WGygoR4kwmE88//zxbtmzhjz/+4LXXXmPLli1leu26detYsGCBn0tYeVpr7H78LAr1QDEfuMYz+qk3kKG1PhjsQp2U73Lxj/RVCBHSGjZsSLdu7gGTMTExtGvXjv37T7/f/PTTT+nYsSNdunShf//+5Ofn88gjjzB37lySkpKYO3cux48f5/zzz6dDhw7ccMMNBX0COTk5DBs2jC5dutCxY0fmzp0b0J8xX2u/pjkPatOTUmoOMACoo5RKBaYAYQBa65nAAuACYCeQC1wXnJIWb0deHk0jIjBJVlkhQt6ePXtYu3YtvXr1Om3f1KlT+eGHH2jcuDHp6emEh4czdepUUlJSePXVVwGYOHEi/fr145FHHuG7777j7bffBuD777+nUaNGBetdZJSS4cHX8l3+6p1wC/aop9Fa64Za6zCtdROt9dta65meIIFntNN4rXUrrXUnrXX583L4mc3lkhFQQlQB2dnZXHrppcyYMYPY2NjT9vft25cxY8bwf//3f8Um/1u6dGnBkqrDhg0jISEBgE6dOvHTTz9x//33s2zZMuLi4vz3gxTi1Bqnn5vA5TbYB3bl5WH3c0QXQlSc3W7n0ksv5aqrruK///1vkcfMnDmTxx9/nH379tG9e/dT0pCXpnXr1qxZs4ZOnToxefJkpk6d6quil8rftQmQQOETdq3ZJbO1hQhJWmuuv/562rVrx1133VXscbt27aJXr15MnTqVunXrsm/fPmJiYsjKyio4pn///nz88ccALFy4kLS0NAAOHDiAxWLh6quv5t5772XNmjX+/aE8tNbkB2BATbUZHhtsu61WmkdGEmk0BrsoQoS00oaz+try5cv54IMP6NSpE0lJSQA8+eSTpw15vffee9mxYwdaawYNGkSXLl1o1qwZ06dPJykpiUmTJjFlyhRGjx5Nhw4d6NOnT8FQ/I0bN3LvvfdiMBgICwvjjTcCk0TC7udO7JMkzbiXsqQZL0m98HB6xsSg/DSWWYiqSNKM+0+2w4HD6zM8xmTCWIbPn/KmGZcahQ8dyc/nb6s1YMulOrXmmN3O0fx8cl2ugk6tSIOBhLAwEkwm4k0mDBK4hKh2nFqfEiT8SQKFj23NzaVOWBixJv+9tcc8AemI3V7shL+TqdAjDAaaR0ZyRmQkETKEV4hqIxCd2CdJoPAxl9b8mZVF//j4MlUBy3PeVJuN3Xl5ZJUjx5TN5eKv3Fx25OVxRkQEbS0WmfMhAk5rLU2yPlSZmdgV6W6QQOEH2U4nf2Rm0j06utKd2zaXiz1WK3us1krdQbi05m+rlYP5+XSIiqJRRESlyiVEWUVGRnL8+HFq164twcJHHFpXKH2Q1prjx48TGRlZrtdJoPCTE56O8eSYGBLCwsr1WpfWHMrPJ9Vm40h+Pr5shbS6XPyZlcXB/HySoqN9WusRoihNmjQhNTWVo0ePBrso1YbN0ydZWKTBUGqfZGRkJE2aNCnX9SRQ+JHV5WJFZiZnms00DA8vtt/CpTXZTqe7Y9pu57jd7veZlgdsNnKdTnrExMiQXuFXYWFhtGjRItjFqDasTic/p6UVeQM5MD6eaD/0j0qg8DOX1mzPzWV7bi4Wo5FaJhMKd5pypydAZDudPq01lFW6w8GyjAx6xcb6tfNdCOE7e222gH9eyKdDAOU6neSG2GJHVpeLPzIz6RsXR5TULIQIaVrroGSsluEvApsnWNgkX5UQIe2I3Y41CP+nEigE4K7trMzMxCHBQoiQtSdImaql6UkUyHA4WJudTY8iUjALUZhTa3I8fWz5LhcOz0zhaKOR2mFhmKUp06dynU6OeCbSBpoECnGKQ/n57MnLo3mA0pCIqkFrzVG7nRN2O5lOJ1ll6G+zGI2caTbTLCJC5k/4wN4grnsjgUKcZnNuLrX8nIZEVA0Ol4u9Nht/W63lHoiR63SyITubVJuNzlFRxMjfU4W5tGavzRa060sfhTjNyTQk/p7LIUJbntPJ0owMNufkVGq03gm7nV/T04N6R1zV7bfZgjrYRAKFKFK208nmnJxgF0MEidXp5PfMTHJ8NJxbA+uzs9mRm+uT89UkWmt2BnlhtKAGCqXUEKXUX0qpnUqpB4rYP0YpdVQptc7zdUMwyllT/WO1cjRInWcieKxOJyt8GCS8bcvNZWN2dkAW26kuDufnkx3k+VdBCxRKKSPwGjAUaA+MVkq1L+LQuVrrJM/XWwEtpGB9drYMma1BtNakZGX5JUictMdqZZ0EizILdm0Cgluj6Ans1Frv1lrnA58AFwexPKIIeS4XW6S5oMbYZ7OR5nD4/TqpNhtrs7MrlAG1Jjlhtwfk91GaYAaKxsA+r+epnm2FXaqU2qCU+kwp1bS4kymlximlUpRSKZKl0rekCapmsLtcbA3gTcF+m401WVkSLEoQCrUJCP3O7G+A5lrrzsBPwHvFHai1nqW1TtZaJ9etWzdgBawp1mVnY6+CTVB2l4s0u539Nhu78vI4lp8vH0zF+Cs3N6CrpoF7JcaVmZlV8m/L3zIcDg6HyA1aMAc27we8awhNPNsKaK2Pez19C3gmAOUSRbC6XGzKyaFrTEywi1ImOU4nu/Ly2GeznRYYTEpRLzycthaLJEL0yHQ4gpYe4pjdzm+eLMYW+X0U2BRCow6DWaNYDSQqpVoopcKBUcB87wOUUg29ng4HtgawfKKQVJuNg0Gc9FMWVqeTNVlZ/JKWxj9Wa5G1B4fWHLDZ+DU9PSiZOEPRlpycoKS6Pynb6WRZRgbH7fYgliJ07LfZOBFC70XQahRaa4dSagLwA2AEZmutNyulpgIpWuv5wESl1HDAAZwAxgSrvMJtQ04OCSZTyC12dDL98tbcXBxlbFpyas2G7GwO5+fTLTq6xq4lnuFwcDQEPpTyXS5WZGRwptlMG4ul1JXaqiuHy8WWEKpNQJBTeGitFwALCm17xOvxJGBSoMslipfvcrEuO5tesbEhk78n2+FgXXZ2hUeHHM7P5/fMTHrHxhJWA4PFrhDpMD1pZ14eR+x2ukZH18g0Mjvy8oKSSrwkNe+3UM3lZGaSumMHJw4dIu3wYZTBQNPERJq2aUNcnTo+ucZRu52tubm0j4ryyfkqSmvNbquVbbm5le6gTnc4WOEJFhE1KFjkOZ0cCMHmxEyHg6Xp6TSPjKSNxRLwAO7SmnyXC5snI64CFBBuMPi1Xyvb4WB3CDaHSqCoBmy5uaz64Qd+/eIL1ixahKOYZoQGzZsz+JprOO+qq4itVatS19yVl0ecyUTjiIhKnaeicp1O1mZnl9qO67DbOXHoEMcOHCD9yBHOaN+exq1aFXlspsPB8owM+sbF1ZhgsdtqDWrfREk08LfVyv78fNqYzTSLjPR5c5Td5SLD4SDD6STL4SDH5SLH6Swxr1KYUsSbTNQLD6dZRITPmiwdLherQ3S4sKqOsyOTk5N1SkpKuV+XbrezLCPDDyXyD6fDwQ/vv8+cZ58l49gxajVowNkjRtCpXz9qN2xIrQYNcOTns2/7dvZt386q779n4/LlhEVEMPCKK7hm8uRKBQyDUvSNjSU+LMyHP1XpDthsbMjOxl7C3+6hf/7hq9dfZ9HHH4PVign3yI0soGGLFvQ4/3wuGDuWRi1bnvbaeJOJPnFxGEOkac1f7C4XP6WlVZnkjxajkUSzmaaVSFvu1JpjdjtH8/M5ardXOjVGmFK0MptpHhlZqVqP1prVWVmVHg47MD6e6Ao21yml/tRaJxe5TwLFv6pSoFjzyy+8/fDD7Nu+nY59+jDynnvo2KcPxlKqxf9s3cp3s2fz04cfEp2QwM3Tp9N3+PAKlyPCYOCs2NiApJB2ac2mnBx25+SwZtEidm3cSKtOneg2aFDBz5125Aizp0zh7y++4AWtudBgIMLrw2Du8OE8m5PDxmXLMBiN3PTMM5wzcuRpHzwNwsNJjokJmX4Yf9iRm8u2Kjjr3mI00jQigqYREWVaHMnmcnEkP5+DnuDgjzv2SIOBrtHR1AkPr9Drt+XksKOyfUVaMzAhQQJFWVU4UKSns+2HHzg6aJAfSuUb+VYr7z72GN++9RaNWrZkzKOP0mvIEJRSmDIzaTJvHvFr12Lev5/jZ53F9vvvB+CM2bPJ6NKF9O7dAfh782Zeuf12dq5fT58LL2Tiyy9jqeAciXCDgd6xscT5MVhYnU5SsrI4ZrMx5fLL2b5mDbbcXCIsFlp368Zjn37K1uXL+XjcOLbn5HDp//7HmwsWcOz887E2bozLZMJgt3NgxAhymzcn6v33MTz+OKPS0mjz3/9y63PPnfbzt4iMpGN0tN9+pmByac3PaWlVep10BdQOCyPeZCLWZCLaaMSlNXbPV5rdzjG7nawAJtRLNJtpXc4RW3utVtZnZ1f62p3uvZd611yDZeDACr1eAkUZWe+8k/BXXuGPzz/neN++fihZ5ezbvp1nx41jz+bNDL/pJq59+GHCPH0E7adM4Yx338WUm0teo0bkNW7MoSFD2DVxIqasLIZ6mliO9+rFzjvu4MigQTidTr587TU+fOopmrZuzeQPP6R+s2YVKluYUvSKjSXBD81QaXY7KVlZWF0uVv/4I8+OG4fVa/hghMXC6H79uOnHH6kfFsbC776jadeuoDUU8w/bcP58uo4fT65SXGezsalbNx7/4gsiCq3s1zU6miaRkT7/mYIt1WplrQ8+nMTp4k0mukZHl3pn79KazTk5PpvoaMrI4D+HD2Pp06dCry8pUNSMHrsysk6eTE7LlnS//nrMqanBLs4p/li4kLvPO48Thw7x8EcfccPjjxPmXc3VmgMXX8zSn3/m5/XrWb5gAbsmTgTAERPD93/9xaYnnsCSmkqv0aPpN3gwMampXHb77Tw6dy7H9u/nnsGD2bZ6dYXKZ9ea3zMzfb44TarVyorMzILhgrs2bsRWqLnkytxc3vzxR9qGh3Nw+nQOHTnCJ88/z+qffsJZzN3kweHDWbpoEc7ERD5zuRiXksLzN9982vEbK7loT6gKxZE11UW6w8GvGRnszssrNkOuzeXi98zMSgeJ8KNH6XjffRjy8nDExeHq2bNS5yuO1Ci8pNvtrFmzhrPPO4/c5s1Z/u23OC0WP5Sw7LTWfPHqq7w/bRpnJiXx4HvvUbthQ8z79pF0223suOsujvXvX+bzqfx8mnz2GS3ffJPfP/uMfE9erNQdO5h65ZUcP3iQSe+8Q/J551W4zI0iIugcFVXpzr0tubnsLtRuW7hGMQV4FFhVrx5Hf/mFB2+5pchmqeL6blR+Ph0mT6bFO+8wFAi/6SZuePzxU45J8HRuV5cJYMftdlZUkb64qi7OZKJeWJh7aWGjkeMOB/ttNp/0lUQePEjvyy7Dsm8fy+fPJyMpyW+d2VKjKCSnVSvWzJxJ7KZNtJ8yJahlsefn8/Ltt/Pe1Kn0HT6cJ7/+mtoNGxK3di39Bw0ift06wtLSynVOHR7Oviuv5NfFi8mvWxflcND62Wc5o3Fjnvv+e5q1acOT117LHwsXVrjcJ9NjVDTdh9XpZGVm5mlBAqDboEG07taNCIuFMGAA8El0NAf//JOV69ezfc0arDk5aK2x5uSwfc0a1ixaVOy1dHg4m55+mlUffkjYjTcy/803WTB79inHpDkcbK+Cnb7FKep9Ff6R4XCwIy+PlZmZ/JSWxhrPyKbKBgnLnj30uegiIg8cYOWcOWQkJfmmwMWQQFGEI+efz/oXX2TnhAlBK0NudjbTrrqKRXPmMPLuu7ln1iwizGZqL19On0suwR4by9JFizh4cQWX8PDc7SekpND62WfpdcUV1DKZmPbFF7Ts1Imnx45l+fz5pZykeHkuFylZWSzPyCCtjOkhtNbstVpZkp5ebEoJo9HIIx99RJu2bbEDL119NZE7dmCIjCyyWcqWm8vujRtLvrBSHB48mLHTpjGyd2/yJ0/mwK5dpxyyMy+P9BBIc1FZOU4nh0IkI6momOi//qLvhRcSlpnJ719+GZD+VAkUxdh31VXknXEGaE34kSMBvfaJQ4d4cPhwNixbxsSXXuKqBx7AYDAQvW0bvUaOJK9JE1Z88w05xUwcK9e1evfmz1mzSFizhj4XX0xtq5Wpn31GYteuPDtuHCu++aZy5/dkBv3Nk4CvqHTSDpeLgzYbv2dmsr6U+REup5OwgQN5fc0abn3sMca++CJGT1W7VadORBRqKoywWGjZqVOZymo0GnmiTh3esNv5adw4XF5l1cDmalCrCFaGWOE72mTCHh/Piq+/9ntN4iTpo/BS1DyKTvfcQ51ly1j28884ApBie9/27Tw2ejSZx45x/9tv0/3cc//dqTWtXnmFfVddRX7t2j69bt3Fi0keMwZbvXqsmD+ftJgYplxxBTvWruXB996jx/nn++Q6BqWIMhiI8HzZXC5OOBxlqoprrck/91wu27CB+b17owoFMafTedrQ2c7Jybzy1Vc4lMLqmXVbkvDjx+nVqxeHMjKYdPXVWJo0OWWuRreYmKDNRq8sm8vFoio0wU540Zp6ixZxZNAg90g+pxOK6HeTCXfl4MtAUWvFCvpccgkHL7yQP996q9jhlr6wacUKnrz2WkxhYTz88cckeu4W4tatwx4fT27z5n67NkB8SgrJN9xAyttvk969OzmZmTx86aXs2bKFhz/6iK4DBvj1+qVJv/xy/rdkCQvbt8e+eDGqiM5yk9YcXraMfzZvpnvXrgwdOvSUjuwsh4NUm41Um63YxGvxv/5K38suYzZwIxAZFVXQKR4dFsbAhIQqOWt7S06OTxMAulwuDu3Zw95t24iwWKjfrBl1mzQ5dTSeqDRTVhYdH3iApvPmsfqddzh04YXFHiuBohx8PTO71csv037aNDY++SR7brzRF0U8zdIvv2TGhAk0OOMMHpkzhwZnnAFA9Pbt9L3oIrLPPJPl337r10AFYLDZcHnumJXDQWZWFg9dcgkHdu/m0Tlz6Bik+SV777iDWz76iGVNm5KxejWqiLupphERdCjjaCuX1uzIy2NHbu5puY5W//gjUddey/0OBxcAC3EHi3tnzaLH+efTxmKhdZBHw5WXL2sTG5cvZ96LL7I9JYW8QumwDQYD7Xv3ZvA113DWsGGEV8M5KIFUZ+lSutx+O+YDB9h+991sv+eegv7FovgrUEhSwDLYNWECtVatosOUKaR37Up6cpHvZYW4XC7mPv88c555hva9e/PQ++8Tk5AAeIa/XXEFLqORta++6vcgARQEiZavv06DhQv5Y948pn76KQ+NGMHUK6/k0XnzaN+rl9/L4W3xp5/y3UcfkdiwIdZffz0tSCigW0wMjcrRJGRQijYWC/XCwlibnX1Kk9SujRv51OHgIHByvNTJTvEe55/Pzrw8mkVEhNyaHCXZlZdX6SCxe+NG3ps2jbWLF1OnUSPOGTWKlu3bk9iiBTkGA4f27WP/zp389vXXPH/zzcQkJDDshhu49LbbTpvIKEp35owZtHviCbJbteK3777z6edOeUmgKAuDgXWvvkrP0aMx+jAlc252NjPGj+ePBQs4Z+RIbn3uuYI7MFNWFj1HjyYsPZ3l8+eT26KFz65bFtaGDam1ciXdbrmFlLffZtrnn/PQiBE8NmoUj336KW0D9Ee7/vPPeWnCBDqefTZZH39c5B1q5+jocgUJbwlhYfSPi2NFZiYZnvUsWnXqhDEqipc8d8sGIMxsLugUd2rNLquVDkFOs15WNper0p3Y37/3HjPvu48u0dE8cf31tJ4yhQizmQuaNsVotWKPjiazQwcyOnXiwIsvssTl4rvZs/nk2Wf5Ze5cxj35JD0HD/bRTxQcWmtSd+zgr5QU/vrzT/Zu20Z2ejpZaWnkW61ExcURHR9PfJ06NGvXjpadOnFmly40SUzE5XIVm5/MW/jRo7giI3HExJDVti27br6ZvyZNCvp8Lml68lJqUkDvlBAlpIcoi73btvH0DTewf+dOxj72GBeNG3dKAroOkyfT/O23WfXxxxytYO6Wymrx5pt0nDyZv6+/nk1PPcXxQ4eYNHw4mceP8+i8eX4PFuu+/JJLxo1jTd26RKxahaWIvEvtLBbO9ME/kc3l4reMDHKdzlM6xdvm5DAPuK1FC27+/feCf26jUpybkEB4FUhHvjknp8JzJ7TWfDR9OjteeIF34+LonJFBRqdOLP3lF8DdLKuNRix79xK3cSOxmzeza/x4tt93H8rhYPvPPzPj8cfZ99df9B42jAkvvFDpFPeBlrpzJ0s+/ZQln33Gkb17AYiKjaVFx47E1q5NdHw8EWYzORkZZKWlceLQIfb+9Rd2z01lQv36oDXZGRnYbbZT+ryMRiPG3Fzq//ADjb/8krqLFrFz4sSCHG3lJX0U5eDX7LFa0/q55wg/epRNTz9d7mDhcrmY/+abfPDEE1hiYrhn1iy6nH32accZs7OptXp10ILESe2nTKHV66+zZcoUdk2YwNH9+3loxAjSjxxh0nvv+a2De/V339HnuusYpDXL3n+fnKFDTzvG10n7sh0OfsvIwK41TqeTNYsWcfCPP3jqtdfYFR7Owd27MXrlsko0m2kb4rWKPKeTX9LTKzTBy+lw8Prtt9Nz3jweVQp7/frsvvVW9l9yCbYGDYp8jcrPx5CfjzM6mgbffkuXO+9ky/3381J2Nh8+8wzxdety1xtv0LGC+YgCRWvNul9/5bMZM9i4fDkGg4Fu/foxondvOjVqhB49GoPBQOLzz1Nr5Up3315kJE6zmbxGjdgwdSqpO3Zg+uAD9v7yC9t37SIPMAPZwKLwcO549VXu++ADaq1ciTE/H2v9+hwYMYJ/rr2W7MTECpW7WvZRKKWGAC/hXjP7La319EL7I4D3ge7AcWCk1npPoMt5CqUw5ubS4p13yE5MLFfn9v5du3j9nnvY+Ntv9BwyhAkvvEC8J4UGAFpzxvvvk3rppTijo4MeJAC2TJlCxJEjODwfiHUbN+bpb79lysiRTLvySu6eObNSacqLsuKbb2h6ww0M0ZpVjz9eZJCIN5l8vsJetMlEj9hYfs/IwGg0uocEn38+S7Oz+e877/DeffdR68UXC47fY7XSymwO6eVTK7P637tTpxI5bx6PA6kjRrDp6aexe/rPiqPDw3F6Rj1lt2lDZseOJE2axIxOnTj/hRe4+4UXmHzJJYy6916uuOsuDCH43v25aBFznnmG7WvWcGa9erw1fDjn2Gw0XrmS8KVLsUdH8/2VVwIQfuIEYRkZuMLDCUtLI/LgQYx5eRhNJs5o146+69ZRq9DkzdXAD/n5PH/zzXRp0oSG/fqhr7mGvCFDihzyGgqCVqNQShmB7cB5QCru92+01nqL1zG3Ap211jcrpUYBl2itR5Z2br+vR+F00uO666j/ww+s+ugjjnjPdShC2pEjfPLcc/zw/vtEREZy/eOPc95VV5221kHrp5+mzXPPsenxx/n7ppvKXX6/8WpmM+bk4IyKIjsjg2lXXslfKSlcO2UKI265pdJrN2it+XTGDGo9+SQvAVtvuIGdTz112nEmpegfH++3JSkLN9Vop5MWLVtS22pl1Y4dGGNjC/b5qunLHzIdDn5NT6/Qa3/55BNm3HYbF15/PQ8MH86JitYAtKbh11/T4ZFHiDhyhHUPPcTdW7ey5NNP6T5oEHe98UbB4I1g2/vXX8x+5BHW/PIL9Zo25bI77mDizp20eeMNrPXrc+Scczh29tlktW1LZseOZWpNCDtxgq0LFvDRgw9izMsjF0gDjkZG0v3cc9m9cSOH//kHg8FA2x496DF4MEn/+Q/NO3QodW2ZogSt6UkplQKsBzZ6vjZorY9WqCSnnvcs4FGt9WDP80kAWuunvI75wXPM70opE3AIqKtLKXRFA8Xlo0djT0igZadOtOzYkSaJiQWzfgszZmfT96KLiN61i1UffMCx//znlP1aa3Zv2MDPc+bwyyefkG+zMfiaaxh5990k1Kt32vlavfIK7adOZe+VV7L+xRdLHAIXLLWXL6f7DTew6sMPSe/eHVtuLi+OH8+Kb7+l97BhTHzpJaLj4ip07nyrlVfuvJNfP/uMu/v147b69dnw2mtF3mElRUfT1I/DLh0uF4vT00+ZZ3H85ZcZM20aH48eTczLLxdsDzcYODdE51X8npHBsQqkHdm3eDEdR47kjc6duWLhQkw+SB1vysyk8z33cGD4cA4OG8YP77/PrEmTqNWwIZNmz6ZVly6VvkZF5WRm8vHTT/PdW29xSUQEzyckkPrMM5wYPJiIQ4eIPHyYjM6dK9wnWdRE0JN9FAaDgV0bNrDq++9Z+f33/L1pEwDm6GjaJifTolMnGrZoQaOWLYmvW5dIi4UIiwXtcpGTkUFOVhYnDh3i4N9/c2D3buopxXuFcpWVVWUDRSOgs+crGRgGHNNan1Gh0vx73suAIVrrGzzP/wf00lpP8Dpmk+eYVM/zXZ5jjhVxvnHAOIBmzZp1/+eff8pVHrvdTp9+/diwYQP5nhEi4ZGRNG/fnpadO9O8XTsatWpF4zPPpHbDhhgMBsKPHqX3yJHsnDCB/ZdcQtrhw2xfs4bta9aQ8vPP7Nm8mbCICPpcdBGj7rmn6LWaXS7aPvUUiTNmsH/ECNbMnBmy1c+Iw4cLcswsnz+f7DZt0Foz/803efexx6jbpAkTZ8wo91yLTStWMPP++8nbto3BDz7I5XfcUWztpFFEBN0DMEP+oM1GSlZWwXOtNXN79mSp0cgry5efcrfXKSqK5iE2/PNIfj4rMzPL/brsnTvpcPbZJDqdLJ89m7wSJndVRtMPP2Sr08nEF14g4/hxbnziCQZfc01AVxTUWvPbV1/x9sMP0+3wYV6pU4e2x46R07w5G595xqdNvyf7vHZv3EjLEkY9HTtwgM1//MHWP/5gy6pVpO7YgaOMubmi4uLo0qkTvy1dWqH30aed2UqpdsBlWutp5S7JqefxaaDwVpmmpyXHj7N/5052bdzI7o0b2b1hA7s3biTH65/OYDBgiY0lKi6OcIOBzKwscjIzqZufz0HAaDKR2LUrAy6/nP6XXEJ0fHyx1ww/doz+Awdy5Lzz2Pj00+gArz9dXpa//6bvhReiXC5+/+wzsjp0AGDrqlU8d9NNHE1NJfncc/nf5Mm08OwrzvGDB3n/8cdZPG8eF9SrxxeZmWx6+WUOXHJJkcdHGAwMjI8PWJ/AysxMjnj9k/721Vc8c+ONTJo1i7O8yhhtNDIgPj5klk3VWvNrenq5V3bTTicJnTvT+8gRvnvmGYzXXeeX8im7nf7nnEP0zp2kPPAAN//2G2uXLKH/f//Lrc8/X+ToNl9L3bmTWZMmsW7JEubWqsUVJ06Q26QJO+66i32jRoXM/6HT6eTY/v0c2L2brLQ0rDk52HJzMRiN7s+g2Fji6tShUcuWxCQkBLXp6Qyt9T+Ftn2gtf5fhUrz7zlCrumpuD4KrTUnDh9m/86d7N+5kxOHDpGdnk52RgYup5OouDi6Z2Ux5auv+GvAAHa/8AKqUaMSrxW7aRM5LVrgjIoi/MgR97oQIfJBU5qoXbs467//xZiXx7IffyxILWLLy+Pbt97is5deIjczkw5nnUXn/v3pcvbZ1G7YkHyrFVtuLltXrWL5N9+w5Y8/MJpMPDpsGPf/8APWBg1Y/u235BfRNAfQIyaGBgHMs5TrdLLYa8SQ0+lkQefOPJKWxsZNm3B4DfPsFRtLvRBJXbHfZmONV22orMLGjWPIl1/y4eDBxH34oR9K9i9TZiZdb7mFBj/+yN5Ro3i0SRPee+EFGjRvzu2vvEI7Py3AY8vNZd6LL7Lw1VfBbObKBx/kpshIzMePs/vmmwsmnFZVwQwUK4BmwN+4+yjSgaFa664VKs2/5zXh7sweBOzH3Zl9pdZ6s9cx44FOXp3Z/9VaX1Hauf3emV0Eg9VK4gsvcOYrr2BPSGDr5MkcuuAC7IVqE1G7dtFm+nQaf/UVWx55hF233Vah6wWb+Z9/OOO999g2efJp/SnZ6enMnzWL1T/8wO6NG4tc5atZ27b0uegiboyOZsC0aWS1bcvKTz7BVr9+kdcLVJNTYZuys/nba7Lapmee4f5nn2XFxReT9tZbBdvrhoXRu4L9M76ktWZxenqpyQ8Lyzx2jKguXTBaLDi3bDllGLDfuFy0fvZZ2jz3HOlJScx86CFevOsujqamMvymm7h60qTTsgFXlNPp5NfPPuOj6dNJTk3lrago/rn3Xo6NH++T84eKoM6jUO46dSugE1AL+OFkc1BlKKUuAGbgHh47W2v9hFJqKpCitZ6vlIoEPgC6AieAUVrr3aWdNxiB4qTYTZvocuedxK9bh7VuXX7avBmUoustt5CQkoJl716ckZHsvukmdo0fjyMEPlwqy7J7N42//JIdd9xxWv9K5okTbFq+nJzMTMLNZsIjImiSmEjT1q2J2byZAQMGcKxfP1a//36x2XnDDQYGxMcTEYQOfqvTySKvWoU9Px995pkMtdlYtmXLKVl8B8THE1PBf1Jf2We1sq4Ca2E/f8stLP/qK2YsXEizAKWuPqnejz8SvXMnu2+9ldzsbN6bOpWF77xDg+bNueqBB+g3YkSFRgCBO3D++fPPvP/EE7g2b+bNuDiGZGSQlZjIxmefDchaDoEkE+7KIZiBAgCnk1qrVxNx+HDBwkKd77oLU3Y2WW3a8M///lds80pVlPj887SdPp1jffuy4fnnS10nI/LgQawNGwLQ8OuvOTxkSIlV/mCn9i48XHblY4/x+Kuv8sdVV3F0xoyC7c0iI+kSgPb14rg8tYlyrfGtNQnXX89133xDj/vuY/S99/qvgGVQ95dfqP/jj3x67rm8Pm0a/2zZQtM2bRh9zz30GjqUsDL+HeRlZ7tzhL31Fvu2b+fu2rV5MjMTg8nEzjvuYOeECegQaSr0JQkU5RD0QFEDNfnkEzrdfz+m3FwOn3ceOydMOGXsvSkjg7qLF9Ng4UIaffMNy3780T0WvRT1w8Pp6TVvIRgKZ17NzcoiIjGRcw0Glu7YgdMz8c/gSesRjJoPwF6rlfXlrE00eecdut53Hw/UqUPP9euDniI88YUXaPvUU2S3aMHal17im0OH+PiZZ9i/cydRsbH0HDKEXkOH0qxNG+o2aUKE2YzWmuz0dI4dOMCm5ctZ9+uvbPztN+y5ubTr1Inzbr6Z/8bH0/KTT9gydSp5TZoE9Wf0p2o5M1tUH6mjRnF04EDOeO89mr/7LrVXrOBEnz6Y//mHsy69FHNqKgank/z4ePaMHUtOGdbWMClF5xBIkRFhMNA8MrJgLQdLTAzfX345T3/6KTemp1PXU0aX1uyxWmkThAl4Lq3Lva63ee9e2k2ezE9A+IwZQQ8SADvuuosTvXqRNHEiZw8fzpkXXMCQWbP49dAhln/zDSsXLmTxvHkFx0fHx5OXnY3Tk9AR4IxmzZjWrRtjd+/m+MCBbLviCtKAP320+FZNJDUKL1Kj8A2DzYay23FGR2Pet4+2TzxBXpMmHD7vPNK7d0eX8Y6nS3Q0zUJkPYN8l4ufvWoVR/btY1xyMiPGj2fMI48UHBfmqVWYAlyrKPdIJ5eLHhddhGXVKkb27s0t8+eHzPBecE9obTlzJq3eeINNTz5J6siRhB8/jjM3l00HD3J4716O7N3LiUOHsHiGiHbMzmbwX39x5i+/EJaZSUaHDmx78EGOBClAhBsMtIyMJN5kIspoJNJgINPh4ITDwTG7nSP5+aethVJZUqMQVYYrIgI8bcl5TZuydubMcp+jTlhYyAQJcP/Te9cq6jVtytkXXMBFM2fSuFYt9k9wT/+xa80/NhutAjwBr7zZYZvOmUODVau4USmGPfNMSAUJAGd0NDvuuYe/b7yxoGmv6Ucf0X7aNM5u3Zr82rVxWCw4o6L405OcM2nCBBr++CMHL7yQ1Cuu4NjZZwclw4FJKVqZzbSMjDzthiE+LIz4sDBams1YnU7+sdn4x2rFVsxqi6FCAoUIOUal6BzETuHitIyM5G+rtWAE1LDx42n27be0fOEFDt54Y0GH/K68PFpERmII0Idvmt1OulfTS1ms6tCB/1OK/ddey/B27fxUssrzHhV46IILAKi1ahWmrCwijh9HHTpEWFoa9lq12H733WycPh1nEP92GoaH0zEqqkyLWkUajbSxWEg0m9ljtbIjL4/8EA0YEiiqMaNSxJlMmJTCCNi05kQFcv8EWluLxW8J/yoj0mikaUQE/3jmVbRNTubD1q05Z/t2Gn3yCanXXgu4O7/3Wq0BS+uxu5yLEimHg3dmzGB9VBSzKrjuQTDknHkmuyZOZFcx+wO9uJe3CIOBTlFRNKzA6DyDUrQ0m2kWEcEuq5WdeXkVzvjrL6GXeU5UWqTBQFuLhXMTEugbF0ev2FiSY2PpGxfHeQkJdIyKIjoEP4gBEkwmWoRQk1NhrcxmvOsJ9e67j9VA02efBe/lVL1qHv6U53RysByrLtZdvJg+PXpw9LvvuGT8eOLq1PFj6ao/BbQ0mxkYH1+hIOHNZDDQxmKhf1wc8UGej1OYBIpqppXZzKCEBBItliJXX4s0GmlhNvOf+HjaR0WFVNZTg1J0iY4OufZyb1FG4ynLrvYeNozXa9emzuHDNJo/v2B7rtPJfh8um1ucPVZrmTtEDXl5dLrvPrKOHyezVi2G33yzX8tW3dULD2dAfDwdoqJ8mn8sxmSiX1wcbS0WQuU/QQJFNWFUim4xMbSPiipT27jB0+E2MD4+ZHIUJZrNQZ/ZXBZnejUpGU0mDBMm8BiwvtBQ3q25udj92Obs1LqgGawsznz5ZaL27GFsXh7D77wzIMn3qqN4k4mzYmPpFRtb4RFGpVFKkWixcJYn8WiwBb8EotIiDQb6xcVVaPay2WikV2wsHcoYYPwl1mQ65QM4lMWaTDTwCq7n/u9/TLdY+PDrr085zuZysb2Ca1WXxX6bDXsZm7cse/Zw5iuv8G1CAhsaNWLomDF+K1d1ZTEaSY6J4ez4eOoE6OaqdlhYSDRFSaCo4kxK0Ss2lthK/iG1NJvpGxuLJQh9Fwr3YkTBDFTllegV1KLj4hg0ejS2zz+n+QMPnHLc33l5ZJZzRFJZlac20fydd3AAN6WlMeqeewgP4X6gUGNUirYWi0/6ISrCbDTSNy6OpkFMYyOBogozKEWPmJhKB4mT4j13Lw0C3BR1ptlMXBVocvIWHxZ2SpPdRTfeSA+nk05vv02t338v2K6BjTk5Pr9+ejmHxG6aPJkRjRtDy5YMGj3a5+WpruqGhXFOfDyJFktQb2QMSpEUE0OHqKig9FtIoKjCkqKjfV4FDjMY6BEbS7sAdaRFG420DtE1p0vjXato1KoVa849lyNKceazz55y3Am7nX3lHMJamn/K2FFusNkIO3GCFd99x/e7dzP6vvuKXd5X/MugFB2iougVG1umORGB0tJspndsbMDziUmgqKJamc1+zah6pqcjLdKPf5BVscnJW62wMOp4rdswZMIEntWa+suWEf/nn6ccuyknp9xrRBTH7nKVeURVyzfe4Jzevfn5ySdp2qYN/UaM8EkZqjOzp8+vpdkckiPw6nhGWwWy5i+BogqKMRppG4C78NphYfSPjz/lw9CXWlssJITIkpMV5V2r6NinDz926sQJg4HE55475TiH1qzJyvLJ3Ir9NltBzqmSRO7fT+KLL/JXs2as+ftvrrzvvgqv61BTxHmGpoZ6U2i4p+bfNTqasAAEMwkUVYzBMww2UHfhEQYDvWNjfd48VCcs7JQP2aqqTng4tTzBTinFkIkTmeJysTE2Fgp9mKc7HGwtZ4bXouwpYzNW+ylTUC4XY48do0WHDpx14YWVvnZ1VjcsjD4h1tRUmiaRkQxKSKC1xYLJj58JEiiqmNZms886r8tKKUUbi4XesbE+GdMdbjDQNcQn1pVHG6+Ad9aFFzKvWTMm7ttX5Brou/PyOJyfX+FrHcvPJ6sMTVi1ly2j8ddf8/2AAaTs389VkyZhCIHx+KGqcUQEPWNjA5711xfCPDO6ByUk+K2puOq9KzVYfJDnGtQND+c/cXHUrmRzUbfo6Cp111aaOuHhBe+J0WRixM0389fq1bief57YTZtOO35tVlaF+yt2lnFeRr1Fi8hp1oyb16+ndbdu9PBxqu1aYWE0ioigQXg4dcPCqmw/E0DTiAi6VuG+spPCDQa/BbqgBAqlVC2l1E9KqR2e7wnFHOdUSq3zfM0v6piapFNUVNDvwiONRs6KjaVNBUZFKaBrdDR1Q2QmuC95L1Z07pVX0jA+nv88/zztH330tGPtWrM6MxNHOWdtp9vtHC1jUsetjz7Kg1ddRerBg/xv8mSf/d3Emkz09uQN6x4TQ4/YWHp7coi1s1j8OvjBH87wLF8b7P+rUBes3+oDwCKtdSKwyPO8KHla6yTP1/DAFS/0NImIID5EOn6VUrS2WOgbF1fmCXoGpegeE0OTajrRq3ZYGHU9v5/IqCj6jR3LY3Y7dX/9lbqLF592fJbTybpyLltaltpE5MGDRO/YQW5WFu/OmkXSgAF0Ofvscl2nOB2iougfF1dkoA83GDjT0/zRMSoqIB2sldUiMpLOEiTKJFiB4mLgPc/j94ARQSpHlWBUinYhONcgISyMAfHxtLVYSvxgMHkmBgZjVmsgeXf4Dx83jtlmMwctFtpNnQpF1B4O5uezo4yd21kOBwfL0LfR6f776Tt0KAtfeonM48f530MPlf0HKMbJmmBZhosalKKFJzFli8jIkElqV1ii2UxHyXVVZsEKFPW11gc9jw8B9Ys5LlIplaKU+kMpNaKkEyqlxnmOTTl69Kgvyxp0iWZzyLbpGz3JywYlJJBoNlMrLKyg+SHOZKJzdDTnJSSETOJBf6rlNVs7tnZtzrvxRu7NzSVu0yYaf/FFka/ZlptLahlGMZWlNtFgwQIaLFzIphtvZN5bb9HnootITEoq189QmEEpkitQEwwzGOgYHc2A+Hjqh9jvvq3FQtsQWIu9KvHb8Bml1M9AgyJ2nXKLo7XWSqniBoWfobXer5RqCfyilNqotS5y3RKt9SxgFrjXzK5E0UOK2WCgZRUYRhpmMJzyz+fUOqRSmAdKG7OZI547/0vGj2fcW2/xa1wc5piYYl+zLjubcIOh2GCa7XCUOsHOmJ1NxwceIKNDB6alp2PLy+PqSZMq/oPwb4qYygT5aJOJnrGxHM3PZ0turt/yXpWFAjpGRQVsQanqxG81Cq31uVrrjkV8fQ0cVko1BPB8P1LMOfZ7vu8GlgBd/VXeUNXWYqmSH7hVscy+4J0DKrZWLYaNG8fAQ4dY1axZsa/RQEpWFmlFdFRnOhysyMwsdc2Jtk88QeShQ/w8cSLfvPsu5119NU0SEyv8cyige3S0z2qCdcPD6R8XR5fo6ICnnwD332OP2FgJEhUUrKan+cC1nsfXAl8XPkAplaCUivA8rgP0BbYErIQhIMZo9GuaDuEf3vMqRtx6K5FRUXz69NMkPv88MZs3F/kap9asyMxkU3Y2eZ6hs8ftdpZnZGArbXSU1jjNZnaPG8cTc+YQabFUqjahgO4xMTTw8d+eUopmkZGcEx/PmYVWCvSnSIOBvnFxIdcEVpUEK1BMB85TSu0AzvU8RymVrJR6y3NMOyBFKbUeWAxM11rXqEDR1mKRERlVUHxYWMGHUkxCAheNG8em776j2cyZdLnzzlOWTPXm0pq/rVYWpaWRkpnJH5mZOMqS8kMptj3yCO/268faJUsYfd99FV7i9OTMf38OPDAZDLSLiqJ/fHzBrHZ/OZmGJtRTcoQ6pUNsEW9fSE5O1ikpKeV+XbrdzrKMDD+UqPziTSbOjo8PdjFEBWU4HCxNTwcgOyODW846i7ExMbywezebp05l9y23VP4iWtNx0iQOXnQRh5KTGd+vH2Hh4by0ZAmmCnwAmzwd14Gc56K1Zo/VypbcXJ+vMZ5oNrvn+8jNVpkopf7UWicXta9qzY6pQUJxOKwouziTqaBWER0Xx5iHH+bF3bvZ1LEjbZ96iujt2yt9jaYffUSLt98mfu1avp45k0N79nD9tGkVChIRBgN9ipkj4U/KM5y2f1ycz1LTRHryk7UNgQmq1YUEihBUJywsYEstCv/xnq09cORI2iQnc/mBA9gtFrpff32xTVBlUXv5cjrdfz9Hzz6bxYMGMefZZ+l9wQV0O+eccp8rFDKmxphMnB0XV+m+i2aRkQyIj6+Ws/+DSQJFCJLaRPXgXaswGAzcNH06f6WlMbl/f9a/9BJUcG5MzObN9Pjf/8ht3pyVM2fy4oQJmKOjuaXQgkll0TQiolwz7P3JoBTtoqLoGxdHTDnLczK1SJfoaMKqWBqRqkDe0RDTOIRSdYjKa+01AurMLl0YfM01vDh/PimeD8Lay5YVOWu7JGe8/z6O6GhWzp3Lh7Nns2vDBsY//zwJ9eqV+RwKd+6wpJiYkBvKnODpgG5nsWAu5UM/ymikW0xMsalFhG9IZ7aXYHdmG5RiYHx8SNzdCd9ZmZlZMAkv88QJ7hg4EIPJxIdPPcV5V13FoSFDWPfSS9hr1QLA6XSyZtEidm3cSKtOneg2aBBGoxFDXh4usxmcTiIPHWL9oUPcP2wYAy67jDtefbXM5VFA15iYKjH0WmvNEbud/TYbdq1xaI3C3TzbIDw84Cn3q7OSOrMlUHgJdqBoaTbTQVILVDtpdju/ef1dbUtJ4cHhw0nq35/Z//kPHaZNw1a3LmtmzuRoz55Mufxytq9Zgy03lwiLhXZdu/JRnz6cMWcOSxctwp6QwPGDB7l/2DC01rz8669ExcaWqSwGpegeHe3zORKi6pNRT1VAmFLVYsU3cboEr8yyAG2Tk7l+2jRSFi3i6dxcflu4EGdEBH1GjKD2DTewfc0aDDk5XKw103JyeH/5cto98wwnevZEG41kpaUx5YoryDpxgknvvlvmIAHQww8T6UT1J/W2ENHaYvHJ6nEiNLWPiuJXz7wKgAvGjmVbSgofP/00dRo3JvuXX0h84QW27NqFLTeXRsCXgBVYrTWvDR9Os5kzsebkMHXkSA7s3s2jc+dyZpcuZS5DotlcI5IzCt+TQBECaoWF0aKartMg3GJNJs6IjOQfT6ZYpRTjn3uOE4cO8dJtt5G6fTtXP/QQ+xctIuLXX9mfk0NPYD1giIri3pEjMaem8sKtt7JjzRrue/ttOvfrV+br1woLO2W4rhDlIbewQWZQii4yMahGaGOxYPL6PUdGRfHYvHkMGTOGz195haeuvZZm7drRuls3IqOiSFEKQ1QUiV27cnjvXib068eu9eu547XX6HPhhWW+bphSdJMFekQlSGe2l2B0ZrePiqKV9E3UGLvz8tick3PKNq01C2bP5v8eegiX00nbnj1p2ro1x/fvx+VyceLwYfZu20bSgAGMf/556peQibYoyTVg0ShReSV1ZkvTUxAlmEy0lCanGqV5ZCR7rFZyvGZlK6UYdv31dB04kGVffslv8+fz04cfAhBXty7N2rTh9lde4ZyRI8tdK2gUESFBQlSa1Ci8BLJGEWU00ic2NmRXrhP+k2a383tmJs4S/vcO792LOSqK2Nq1K3ydMKUYmJAQlPUfRNUjw2NDjMVo5CwJEjVWQlgYPWJiMJRQO6jfrFmlggRAh6goCRLCJ+SvKMAiDQbOio3FLEGiRqsbHk5SdLTfzl8nLIym0qwpfET6KAKoUUQE7S0WCRICoCCFxo7cXLIqkUm2MJNSdPFjEBI1jwSKAIg1megYFUVtSfYnCmkcEUHjiAgyHA7222zkuVw4tMbucpHucJS6VnZRusXESL4w4VMSKMopTCmUUji0LnJFLpNSRBgMRBuN1PWkboiWxGWiFHEm02nrQeQ6nezOy2OvzVZix7e3thaLrA0tfC4on2BKqcuBR3Gvi91Ta13kECWl1BDgJcAIvKW1nh6wQnrUCQujpdlMLZMJkydInOTSuuAfWOPu8DFJ56HwEYvRSMfoaNpYLGzKySHVZivx+Ibh4STK7GvhB8G61d0E/Bd4s7gDlFJG4DXgPCAVWK2Umq+13hKIAsabTHSJji4xjbFBqRJHrgjhC2EGA11jYmgQHs6GnBzyC61foYAzIiNlwSvhN0EJFFrrrUBpk4d6Aju11rs9x34CXAz4PVDEeVbLkpWyRChpGBFBrbAwDtpsnHA4SHM4iDMaaWuxSPOm8KtQ/utqDOzzep4K9CruYKXUOGAcQLNypjjwFmM0SpAQISvCYKC52UzzYBdE1Ch+CxRKqZ+BBkXsekhr/bWvr6e1ngXMAvfM7IqcI8Jg4Ky4OEn3LYQQXvwWKLTW51byFPuBpl7Pm3i2+Y3MbxBCiNOF8q3zaiBRKdVCKRUOjALmB7lMQghR4wQlUCilLlFKpQJnAd8ppX7wbG+klFoAoLV2ABOAH4CtwDyt9eZglFcIIWqyYI16+hL3So+Ftx8ALvB6vgBYEMCiCSGEKCSUm56EEEKEAAkUQgghSiSBQgghRIkkUAghhCiRBAohhBAlkkAhhBCiRBIohBBClEgChRBCiBJJoBBCCFEiCRRCCCFKJIFCCCFEiSRQCCGEKJEECiGEECWSQCGEEKJEEiiEEEKUSAKFEEKIEkmgEEIIUSIJFEIIIUoUrDWzL1dKbVZKuZRSySUct0cptVEptU4plRLIMgohhHALyprZwCbgv8CbZTh2oNb6mJ/LI4QQohhBCRRa660ASqlgXF4IIUQ5hHofhQZ+VEr9qZQaF+zCCCFETeS3GoVS6megQRG7HtJaf13G0/TTWu9XStUDflJKbdNaLy3meuOAcQDNmjWrUJmFEEKczm+BQmt9rg/Osd/z/YhS6kugJ1BkoNBazwJmASQnJ+vKXlsIIYRbyDY9KaWilFIxJx8D5+PuBBdCCBFAwRoee4lSKhU4C/hOKfWDZ3sjpdQCz2H1gd+UUuuBVcB3Wuvvg1FeIYSoyYI16ulL4Msith8ALvA83g10CXDRhBBCFBKyTU9CCCFCgwQKIYQQJZJAIYQQokQSKIQQQpRIAoUQQogSSaAQQghRIgkUQgghSiSBQgghRIkkUAghhCiRBAohhBAlkkAhhBCiRBIohBBClEgChRBCiBJJoBBCCFEiCRRCCCFKJIFCCCFEiSRQCCGEKJEECiGEECWSQCGEEKJEQQkUSqlnlVLblFIblFJfKqXiizluiFLqL6XUTqXUAwEuphBCCIJXo/gJ6Ki17gxsByYVPkApZQReA4YC7YHRSqn2AS2lEEKI4AQKrfWPWmuH5+kfQJMiDusJ7NRa79Za5wOfABcHqoxCCCHcTMEuADAWmFvE9sbAPq/nqUCv4k6ilBoHjPM8zVZK/VXB8tQBjlXwtf4k5SofKVf5SLnKpzqW64zidvgtUCilfgYaFLHrIa31155jHgIcwEeVvZ7WehYwq7LnUUqlaK2TK3seX5NylY+Uq3ykXOVT08rlt0ChtT63pP1KqTHAhcAgrbUu4pD9QFOv500824QQQgRQsEY9DQHuA4ZrrXOLOWw1kKiUaqGUCgdGAfMDVUYhhBBuwRr19CoQA/yklFqnlJoJoJRqpJRaAODp7J4A/ABsBeZprTcHoGyVbr7yEylX+Ui5ykfKVT41qlyq6FYfIYQQwk1mZgshhCiRBAohhBAlqpGBQil1uVJqs1LKpZQqdihZcSlEPB3sKz3b53o6231RrlpKqZ+UUjs83xOKOGagp1/n5JdVKTXCs+9dpdTfXvuSAlUuz3FOr2vP99oezPcrSSn1u+f3vUEpNdJrn0/fr9JSziilIjw//07P+9Hca98kz/a/lFKDK1OOCpTrLqXUFs/7s0gpdYbXviJ/pwEq1xil1FGv69/gte9az+99h1Lq2gCX60WvMm1XSqV77fPL+6WUmq2UOqKU2lTMfqWUetlT5g1KqW5e+yr/Xmmta9wX0A5oAywBkos5xgjsAloC4cB6oL1n3zxglOfxTOAWH5XrGeABz+MHgKdLOb4WcAKweJ6/C1zmh/erTOUCsovZHrT3C2gNJHoeNwIOAvG+fr9K+nvxOuZWYKbn8Shgrudxe8/xEUALz3mMASzXQK+/oVtOlquk32mAyjUGeLWI19YCdnu+J3geJwSqXIWOvw2YHYD3qz/QDdhUzP4LgIWAAnoDK335XtXIGoXWeqvWurSZ20WmEFFKKeAc4DPPce8BI3xUtIs95yvreS8DFurihxj7SnnLVSDY75fWervWeofn8QHgCFDXR9f3VpaUM97l/QwY5Hl/LgY+0VrbtNZ/Azs95wtIubTWi73+hopLqeNrlUnRMxj4SWt9Qmudhjt33JAglWs0MMdH1y6W1nop7pvC4lwMvK/d/gDilVIN8dF7VSMDRRkVlUKkMVAbSNf/5qo6ud0X6mutD3oeHwLql3L8KE7/I33CU/V8USkVEeByRSqlUpRSf5xsDiOE3i+lVE/cd4m7vDb76v0q7u+lyGM870cG7venLK/1Z7m8XY/7zvSkon6ngSzXpZ7fz2dKqZMTcEPi/fI00bUAfvHa7K/3qzTFldsn71Uo5HryC1WGFCLBUFK5vJ9orbVSqtixy567hU6455mcNAn3B2Y47vHU9wNTA1iuM7TW+5VSLYFflFIbcX8YVpiP368PgGu11i7P5gq/X9WRUupqIBn4j9fm036nWutdRZ/B574B5mitbUqpm3DXxs4J0LXLYhTwmdba6bUtmO+X31TbQKFLSSFSBsWlEDmOu1pn8twVliu1SEnlUkodVko11Fof9HywHSnhVFcAX2qt7V7nPnl3bVNKvQPcE8hyaa33e77vVkotAboCnxPk90spFQt8h/sm4Q+vc1f4/SpCWVLOnDwmVSllAuJw/z35M11Nmc6tlDoXd/D9j9badnJ7Mb9TX3zwlVourfVxr6dv4e6TOvnaAYVeu8QHZSpTubyMAsZ7b/Dj+1Wa4srtk/dKmp6KV2QKEe3uIVqMu38A4FrAVzWU+Z7zleW8p7WNej4sT/YLjACKHCHhj3IppRJONt0opeoAfYEtwX6/PL+7L3G3335WaJ8v36+ypJzxLu9lwC+e92c+MEq5R0W1ABKBVZUoS7nKpZTqCryJO6XOEa/tRf5OA1iuhl5Ph+PO0ADuWvT5nvIlAOdzas3ar+XylK0t7s7h3722+fP9Ks184BrP6KfeQIbnRsg375U/euhD/Qu4BHdbnQ04DPzg2d4IWOB13AW4F1bahftu9OT2lrj/kXcCnwIRPipXbWARsAP4Gajl2Z4MvOV1XHPcdwqGQq//BdiI+wPvQyA6UOUC+niuvd7z/fpQeL+AqwE7sM7rK8kf71dRfy+4m7KGex5Hen7+nZ73o6XXax/yvO4vYKiP/95LK9fPnv+Dk+/P/NJ+pwEq11PAZs/1FwNtvV471vM+7gSuC2S5PM8fBaYXep3f3i/cN4UHPX/Lqbj7km4GbvbsV7gXetvluXay12sr/V5JCg8hhBAlkqYnIYQQJZJAIYQQokQSKIQQQpRIAoUQQogSSaAQQghRIgkUQgghSiSBQgghRIkkUAgRAEqpxUqp8zyPH1dKvRLsMglRVtU215MQIWYKMFUpVQ93/p/hQS6PEGUmM7OFCBCl1K9ANDBAa50V7PIIUVbS9CREACilOgENgXwJEqKqkUAhhJ95sqB+hHsVsmyllK9WYxMiICRQCOFHSikL8AVwt9Z6KzANd3+FEFWG9FEIIYQokdQohBBClEgChRBCiBJJoBBCCFEiCRRCCCFKJIFCCCFEiSRQCCGEKJEECiGEECX6f9X0jFPuhLChAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAASuElEQVR4nO3deZCkdX3H8fdHEY1CBNkRcUEXSzSiaEmNxitRwSQcUTBSFJbHQm2yVd63oElFc1RlSeJF4lGrqKtRAfEA7xiOwgt0EcOpccVFFoEdoxiPeKx+80c/WxmGnp2e6Znumd++X1Vb0/08Tz/92dmez/7m9zz9dKoKSVJb7jTuAJKkxWe5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdy0KSa5I8adw5llKSByf5RpKfJHnxIu3zxiSPXIx9qS2Wu5Zckq1JnjJj2clJvrjzflU9tKounu9+VphXAxdV1d5VdcawO0uyL3AAcN3QydQcy10Ckuwxgqe5P3DNfB+0i2yHAVuq6hdDpVKTLHctC9NH5UlOTXJTN33xrSRHJnk/cD/gE0l+muTV3bYPSXJxktu6qZ2nTdvn4Umu6Pbz4SRnJ/n7Gc95apIrgZ8l2SPJaUm+0z3m2iRPn7H9q5JcmeRnSc5Msn+Sz3Tb/0c3mu7397sQeDLwr13+B82R/Q7Z+uz24cDV3fZ3T/LBJB9NstfC/yXUCstdy0qSBwMvBB5VVXsDfwJsrarnAN8DnlpVe1XVPya5C/AJ4N+BewMvAj7QzW3vCXwMeC9wL+BDwNPv8ITwTOBYYJ+q2gF8B/gD4J7A3wD/luSAads/A/gj4EHAU4HPAK8FJuj9PPWdS6+qI4AvAC+sqr2A786WfRfZZjoMuCrJwcCXgG8Bz6iqn/bLoN3LKH4VlQA+nmR6Qe0JfL3Pdr8B7gocmmSqqrbuYp+PAfYCNlTVb4ELk3ySXileSO/1fUb1Ln360SRf7bOPM6rqxp13qurD09adneQ1wKOB87pl/1JVtwIk+QKwvaqu6O5/DDhyF3kHzf76ftn6eDhQwEXAS6rqvF1sq92MI3eNyvFVtc/OP8Dz+21UVVuAl9IruO1Jzkpy31n2eV/gxq4cd7oBWN2tu6luf03rfkV5u2VJntud0XJbktuAhwGrpm1y67Tb/9vn/qBTIrvKvqu8O3Omy/Z04O2zFXsSf8Z3U/7Da9mpqg9W1RPoHYAs4PSdq2Zs+n3goBkFdj/gJuBmYHVXgjsd1O/pdt5Icn/gnfSmhfbr/hO6Gkifxw1rV9nvkK2Pg7uvTwFekWRy+soklyd5B72/j3ZDlruWlW6+/IgkdwV+QW80vHN0eyvwgGmbXwb8HHh1krt058k/FTgL+Aq9KZ4XdgdKj6M3vbIr96BXqFNdllPojY6Xwq6yD+LhwJVVdRWwHvjYzmMDSVbRm8d/bVWtW+zgWhksdy03dwU2AD8AbqFXUq/p1v0D8FfdlMkrq+pX9Arx6G77twHPrapvduv+DFgH3AY8G/gk8MvZnriqrgXeQO8/hlvpHbD80mL/BbvnmjX7gLs4DLiy29fHgY30jmvcjV7xf7CqfrjYubVyxI/Z0+4iyWXAO6rqPePOspSSvBTYVlXnjjuLxseRu5qV5IlJ7tNNy6ylN6L97LhzjcBhwDfGHULj5amQatmDgXPozaVfD5xQVTePN9LSc55d4LSMJDXJaRlJatCymJZZtWpVrVmzZtwxJGlFufzyy39QVRP91i2Lcl+zZg2bN28edwxJWlGS3DDbOqdlJKlBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQcviHapaGmtO+9RA223dcOwSJ5E0ao7cJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lq0JzlnuTdSbYnuXrasnsl+XySb3df9+2WJ8kZSbYkuTLJ4UsZXpLU3yAj9/cCR81YdhpwQVUdAlzQ3Qc4Gjik+7MeePvixJQkzcec5V5VlwA/nLH4OGBTd3sTcPy05e+rnkuBfZIcsEhZJUkDWuic+/5VdXN3+xZg/+72auDGadtt65bdQZL1STYn2Tw1NbXAGJKkfoY+oFpVBdQCHrexqiaranJiYmLYGJKkaRZa7rfunG7pvm7vlt8EHDRtuwO7ZZKkEVpouZ8PrO1urwXOm7b8ud1ZM48Bfjxt+kaSNCJzXvI3yYeAJwGrkmwDXgdsAM5Jsg64ATix2/zTwDHAFuDnwClLkFmSNIc5y72qnjnLqiP7bFvAC4YNJUkaju9QlaQGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KD5rxwmEZnzWmfGmi7rRuOXeIkklY6R+6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUoKE+iSnJy4A/Bwq4CjgFOAA4C9gPuBx4TlX9asicK9qgn7AkSYtlwSP3JKuBFwOTVfUw4M7AScDpwJuq6oHAj4B1ixFUkjS4Yadl9gB+J8kewN2Bm4EjgHO79ZuA44d8DknSPC243KvqJuCfge/RK/Uf05uGua2qdnSbbQNWDxtSkjQ/w0zL7AscBxwM3Be4B3DUPB6/PsnmJJunpqYWGkOS1Mcw0zJPAb5bVVNV9Wvgo8DjgX26aRqAA4Gb+j24qjZW1WRVTU5MTAwRQ5I00zDl/j3gMUnuniTAkcC1wEXACd02a4HzhosoSZqvYebcL6N34PTr9E6DvBOwETgVeHmSLfROhzxzEXJKkuZhqPPcq+p1wOtmLL4eePQw+5UkDcd3qEpSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQUNdOExtmM8HeG/dcOwSJpG0WBy5S1KDLHdJapDlLkkNcs59BZrPHPm4ntu5eWm8HLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1aKhyT7JPknOTfDPJdUkem+ReST6f5Nvd130XK6wkaTDDjtzfAny2qn4PeARwHXAacEFVHQJc0N2XJI3Qgss9yT2BPwTOBKiqX1XVbcBxwKZus03A8cNFlCTN1zAj94OBKeA9Sa5I8q4k9wD2r6qbu21uAfYfNqQkaX6GKfc9gMOBt1fVI4GfMWMKpqoKqH4PTrI+yeYkm6empoaIIUmaaZhy3wZsq6rLuvvn0iv7W5McANB93d7vwVW1saomq2pyYmJiiBiSpJkWXO5VdQtwY5IHd4uOBK4FzgfWdsvWAucNlVCSNG/DfkD2i4APJNkTuB44hd5/GOckWQfcAJw45HMsW+P8oGpJ2pWhyr2qvgFM9ll15DD7lSQNx3eoSlKDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAYNe5671Neg7wHYuuHYJU4i7Z4cuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ3y2jJ9+NmoklY6R+6S1CDLXZIaZLlLUoOcc9eK4PXhpflx5C5JDbLcJalBlrskNchyl6QGWe6S1CDPltFY+W5gaWk4cpekBg1d7knunOSKJJ/s7h+c5LIkW5KcnWTP4WNKkuZjMUbuLwGum3b/dOBNVfVA4EfAukV4DknSPAxV7kkOBI4F3tXdD3AEcG63ySbg+GGeQ5I0f8OO3N8MvBr4bXd/P+C2qtrR3d8GrB7yOSRJ87Tgck/yp8D2qrp8gY9fn2Rzks1TU1MLjSFJ6mOYkfvjgacl2QqcRW865i3APkl2nmJ5IHBTvwdX1caqmqyqyYmJiSFiSJJmWvB57lX1GuA1AEmeBLyyqp6V5MPACfQKfy1w3vAxpcF49UipZynOcz8VeHmSLfTm4M9cgueQJO3CorxDtaouBi7ubl8PPHox9itJWhjfoSpJDbLcJalBlrskNchyl6QGWe6S1CDLXZIatNt8WIcfCiFpd+LIXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWrQbvMZqtJ0g36m7tYNxy5xEmlpOHKXpAY5cpd2YdARPgw+yve3Bo2CI3dJapDlLkkNWnC5JzkoyUVJrk1yTZKXdMvvleTzSb7dfd138eJKkgYxzJz7DuAVVfX1JHsDlyf5PHAycEFVbUhyGnAacOrwUSX14xy++lnwyL2qbq6qr3e3fwJcB6wGjgM2dZttAo4fMqMkaZ4WZc49yRrgkcBlwP5VdXO36hZg/1kesz7J5iSbp6amFiOGJKkzdLkn2Qv4CPDSqvqf6euqqoDq97iq2lhVk1U1OTExMWwMSdI0Q5V7krvQK/YPVNVHu8W3JjmgW38AsH24iJKk+RrmbJkAZwLXVdUbp606H1jb3V4LnLfweJKkhRjmbJnHA88BrkryjW7Za4ENwDlJ1gE3ACcOlXAO83kHobSUfC1qOVlwuVfVF4HMsvrIhe5XkjQ836EqSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1aJgLh0nSsrY7fwShI3dJapAjd2mZ8hLCGoYjd0lqkCN3ScvG7jxHvtgcuUtSgyx3SWqQ5S5JDXLOXdKK45lEc3PkLkkNstwlqUGWuyQ1yDl3aTex2PPU8znX3Dny0XPkLkkNcuQuSQOaz28g434XrSN3SWqQI3dJWgLjvk6OI3dJapAjd0kL0tIZMC39XXZakpF7kqOSfCvJliSnLcVzSJJmt+jlnuTOwFuBo4FDgWcmOXSxn0eSNLulGLk/GthSVddX1a+As4DjluB5JEmzWIo599XAjdPubwN+f+ZGSdYD67u7P03yrRmbrAJ+sAT5lpq5R2clZgZzj9qyzp3TZ101SO77z7ZibAdUq2ojsHG29Uk2V9XkCCMtCnOPzkrMDOYetd0191JMy9wEHDTt/oHdMknSiCxFuX8NOCTJwUn2BE4Czl+C55EkzWLRp2WqakeSFwKfA+4MvLuqrlnArmadslnmzD06KzEzmHvUdsvcqarFCiJJWia8/IAkNchyl6QGjb3c57pUQZK7Jjm7W39ZkjVjiDkz01yZX57k2iRXJrkgyaznoo7SoJeFSPKMJJVkWZw+NkjuJCd23/Nrknxw1Bn7GeB1cr8kFyW5onutHDOOnDMyvTvJ9iRXz7I+Sc7o/k5XJjl81Bn7GSD3s7q8VyX5cpJHjDpjP3Plnrbdo5LsSHLCwDuvqrH9oXfA9TvAA4A9gf8EDp2xzfOBd3S3TwLOXgGZnwzcvbv9vHFnHjR3t93ewCXApcDkSsgNHAJcAezb3b/3Csm9EXhed/tQYOsyyP2HwOHA1bOsPwb4DBDgMcBl4848YO7HTXt9HL1Sck97LV0IfBo4YdB9j3vkPsilCo4DNnW3zwWOTJIRZpxpzsxVdVFV/by7eym9c/3HbdDLQvwdcDrwi1GG24VBcv8F8Naq+hFAVW0fccZ+BsldwO92t+8JfH+E+fqqqkuAH+5ik+OA91XPpcA+SQ4YTbrZzZW7qr688/XB8vmZHOT7DfAi4CPAvF7X4y73fpcqWD3bNlW1A/gxsN9I0vU3SObp1tEb6YzbnLm7X7EPqqrldP3TQb7fDwIelORLSS5NctTI0s1ukNyvB56dZBu9UdmLRhNtKPN9/S9Hy+Vnck5JVgNPB94+38d6PfcllOTZwCTwxHFnmUuSOwFvBE4ec5SF2IPe1MyT6I3ILklyWFXdNs5QA3gm8N6qekOSxwLvT/KwqvrtuIO1KsmT6ZX7E8adZUBvBk6tqt/Od8Ji3OU+yKUKdm6zLcke9H59/e/RxOtroMsrJHkK8JfAE6vqlyPKtitz5d4beBhwcfciug9wfpKnVdXmkaW8o0G+39vozaH+Gvhukv+iV/ZfG03EvgbJvQ44CqCqvpLkbvQuFrUcppVms2IvL5Lk4cC7gKOrapwdMh+TwFndz+Qq4JgkO6rq43M+cswHE/YArgcO5v8POj10xjYv4PYHVM9ZAZkfSe9g2iHjzDrf3DO2v5jlcUB1kO/3UcCm7vYqetMG+62A3J8BTu5uP4TenHuWwfd8DbMfmDyW2x9Q/eq48w6Y+37AFuBx4845n9wztnsv8zigOtaRe81yqYIkfwtsrqrzgTPp/bq6hd6Bh5PGl3jgzP8E7AV8uPsf93tV9bSxhWbg3MvOgLk/B/xxkmuB3wCvqjGPzAbM/QrgnUleRu/g6snV/RSPS5IP0ZveWtUdC3gdcBeAqnoHvWMDx9Aryp8Dp4wn6e0NkPuv6R2re1v3M7mjlsGVIgfIvfB9j/m1JElaAuM+W0aStAQsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktSg/wMqGX0LWvH1ZAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.45376086 0.2549733\n"
     ]
    }
   ],
   "source": [
    "u_pred, k_pred = model.predict(x_test, samples, [process_u, process_k])\n",
    "plot1d(x_u_train, u_train, x_test, u_test, u_pred[..., 0], ylabel=\"$u$\", ylim=[-2.0, 2.0], title=\"$LogN(0, 1^2)$\")\n",
    "plt.hist(k_pred, bins=30)\n",
    "plt.title(\"Histogram for $k_r$\")\n",
    "plt.show()\n",
    "print(np.mean(k_pred), np.std(k_pred))\n",
    "\n",
    "# sio.savemat(\n",
    "#     \"LogNormal.mat\",\n",
    "#     {\n",
    "#         \"x_u_train\": x_u_train, \"u_train\": u_train,\n",
    "#         \"x_f_train\": x_f_train, \"f_train\": f_train,\n",
    "#         \"x_test\": x_test, \"u_test\": u_test, \"f_test\": f_test,\n",
    "#         \"u_pred\": u_pred, \"f_pred\": f_pred, \"k_pred\": k_pred,\n",
    "#         \"noise\": noise,\n",
    "#     }\n",
    "# )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b282b6c8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
